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ABSTRACT 


The Navy Navigation Satellite System (NNSS) is over fifteen years old 
and has well established its utility in position determination, from the knowledge 
of the satellite ephemeris and the observations of the Doppler shifts at the 
receiver. 

The satellite ephemeris is available in two forms. The broadcast 
ephemeris which can be received on real time basis is predictive and has 
larger uncertainties than the post fitted precise ephemeris whose availability 
is restrictive. 

This study is an effort to improve station position recovery using broad- 
cast ephemeris in Doppler data reduction. 

Comparison of precise and broadcast ephemerides, treating the former 
as the standard, yields information about the state disturbance that can be 
associated with the broadcast ephemeris. Statistical information about the 
state disturbance has been used with current observational data for improved 
position recovery. 

The rank deficiency problem encountered in the Short Arc Geodetic Ad- 
justment (SAGA) procedure has been analysed and it has been deduced that the 
fundamental rank deficiency is six, scale information being derivable from the 
wavelength of transmission. Coordinate differences between stations coobserv- 
ing a pass are estimable. 

The uncertainty of the broadcast ephemeris, now in the WGS72 system, 
has been assessed. It is conservatively estimated that its p'> jitional uncertainty 
may vary between 19 to 26 m in-track, 15 to 20 m cross-track and 9 to 10 m 
in radial directions depending on the incidence of the epoch of observations in 
the interinjection period. 
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The broadcast ephemeris indicates a radial bias of -5 m, which appears 
to be the consequence of 

(a) a 9.8 m offset of the origin of the coordinate system, with 
respect to the geocentre, in the direction of 90°W longitude, 

(b) a scale correction of -1.4 ppm required to make it compatible 
with a terrestrial system obtainable from a scale corrected 
precise ephemeris. 

Considering the state disturbance as a signal, sample autocovariances 
have been computed for acceleration, velocity, and position. 

Two specific experiments have been conducted. In the first experiment 
in which three stations coobserve 12 passes, the removal of the radial bias, 
resulted in bringing a chord distance in better agreement with ground survey 
though the uncertainties were unchanged. 

In the second experiment, the station positions of the first experiment, 
simulated range rate observations and the autocovariance function of the 
acceleration signal in the Cartesian coordinate system has been used in an 
adaptive filtering procedure to improve the state of the satellite for one pass. 
The position improved by only 0.6 m while the positional uncertainty improved 
by about 4 meters. Much better results are expected with the use of positional 
signal in the polar coordinate system which could also evaluate the in-track, 
out-of-plane and radial biases in the individual passes. The experiment has 
demonstrated the feasibility of this approach. 
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1. INTRODUCTION 


1.1 General Background and Brief Description of Present Study 

The Navy Navigation Satellite System (NNSS) also known as the TRANSIT 
system is over fifteen years old and has well established its utility in position 
determination. 

The operational satellites of this system (currently five) in near polar 
orbit broadcast a pair of signals with a fixed frequency relationship which is 
received at ground stations and examined for Doppler shift due to the relative 
motion between the satellite and the receiver. From the integrated Doppler 
shift position of the receiver antenna can be computed with the knowledge of 
the satellite ephemeris. 

The major factors affecting the accuracy of a receiver position in 
Doppler survey are the following: 

(1) the type of Doppler receiver and the design of the observations 

(2) the accuracy of orbital ephemerides 

(3) the method of data reduction including the corrections for atmospheric 
refraction. 

The satellite ephemeris for the NNSS is available from two sources: 

(a) on a real time basis, broadcast as a message from the satellite, 

(b) in a more precise form, maintained by the Defense Mapping Agency (D M A ) . 

The precise ephemeris is a set of values for earth-fixed positions and 
velocities at ono-minutc intervals computed by fitting 18-hour orbital arcs to 
Doppler data from a world-wide network (Sims, 1972( . The broadcast ephemeris 
injected into satellite memories twice per day and broadcast automatically to 
users is computed by fitting 3 G -hour orbital arcs to Doppler data from the 
four U.S. Naval Astronautics Group (NAG) stations in Maine, Minnesota, 
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California, and Hawaii and extrapolating these arcs 16 hours beyond the time 
of the last data used [Piscane et al. , 1973]. 

The broadcast ephemeris differs from the precise ephemeris in the 
following respects: 

(i) The broadcast ephemeris is available on real time basis and is predictive, 
while the precise ephemeris is a post-fitted ephemeris. 

(ii) The broadcast ephemeris is based on Doppler data from only four stations 
in the U.S. and generated by the NAG (independently of DMA), while the 
precise ephemeris is based on Doppler data from over 20 stations around 
the globe including the four stations tracking the satellites for the broad- 
cast ephemeris. 

(iii) There are variations in the mathematical models , the parameters used, and 
truncation errors between the precise and broadcast ephemerides (e.g. , the 
corrections to semi-major axis and out-of-plane orbit components are broad- 
cast to the nearest ten meters) [Moffett, 1973]. 

(iv) The broadcast ephemeris is available for all satellites, while the precise 
ephemeris is available for only two satellites (which differ from time to 
time) and only after a time lag. 

(v) The precise ephemeris is believed to have uncertainties of two meters 
in each coordinate [Anderle, 1976], while the broadcast ephemeris is 
expected to have an uncertainty of 12 - 28 m. 

There has, therefore, been an ongoing effort to improve station position 
recovery using broadcast ephemeris in Doppler data reduction. The ap- 
proaches tried out so far can be considered to fall into one of the following 
categories: 

(i) Approaches in which the broadcast ephemeris is allowed to adjust 
by assigning suitable a priori variances to the ephemeris in a computation 
in which up to six unknowns (satellite state vector) are also solved for in 
each pass. Fbr example, the short arc procedure [Brown, 1976] falls in 
this category. 
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(ii) Approaches in which the data is examined pass by pass in the 
"Guier Plane" to account for certain biases and to edit the data before subse- 
quent adjustment taking advantage of the fact that both refraction and ephemeris 
errors are correlated between stations which track the same satellite. For 
example, the procedures in [Kouba and Wells, 1976] fall in this category. 

Further approaches appear to be available. If the post-fitted precise 
ephemeris is considered to representthetruestateofthesatellite, the comparison 
between the precise ephemeris and the predicted broadcast ephemeris (gener- 
ated by an independent agency) over a length of time may yield statistical infor- 
mation which when suitably used with observational data, broadcast ephemeris 
and procedures in stochastic filtering theory (where necessary) may give an 
improved state of the satellite and consequently an improved recovery of 
station positions. An effort in this direction is the central theme of this study. 

The study has been carried out in the following sequence: Most geodetic 
problems are intimately connected with reference frames and solutions where 
adequate precautions are not taken would yield values which may not necessar- 
ily refer to the reference frame of interest. One cause for this is the rank 
deficiency encountered in a normal matrix formed without introducing appropri- 
ate constraints and is closely related to the mathematical model used to relate 
the observables with the unknowns for the solution of the geodetic problem. An 
understanding of the rank deficiency encountered helps in ensuring that appro- 
priate constraints are enforced. An analytical study of the rank deficiency 
problem in Doppler survey in the short arc mode which has been used in this 
study has, therefore, first been carried out and described in Chapter 2 along 
with a discussion about the estimable quantities in Doppler survey. 

The next step in the study was to assess the uncertainty of the broad- 
cast ephemeris by comparing it with the precise ephemeris. Studies of this 
nature carried out by Wells [1974] and White et al. [1975] were based on data 
pertaining to the period before the corrp utational procedure of the broadcast 
ephemeris was upgraded in December, 1975. It was, therefore, felt 
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appropriate to carry out a fresh assessment as part of this study. This is 
described in Chapter 3. The data used for this purpose refers pa illy to 
the pre-1975 period and partly to 197(1. 

Considering that the broadcast ephemeris provides the nominal state of 
the satellite which is sought to be improved, comparisons between the precise 
and broadcast ephemerides may be considered to yield information about the 
shite disturbance. Considering the state disturbance as a signal, the pi'oce- 
dures foi obtaining the statistics of the signal by normal sampling methods have 
been described along with the results obtained in Chapter 4. 

Based on the above, two specific experiments have been carried out and 
are described with their results in Chapter 5. The first experiment was 
designed to study the influence of removing the radial bias in the broadcast 
ephemeris noticed in the investigations in Chapters 3 and 4. With JMR 1 
receiver data from three stations, shit ion positions have been obtained from 
12 coobserved passes, both before and after removing the radial bias, and 
the results have been compared. The second experiment was designed to 
judge the feasibility of using the adaptive filtering technique for local orbit 
improvement. With the station positions obtained ia the first experiment held 
fixed, the precise mid broadcast ephemeris shite vectors for one pass of 
satellite 19 and the parameters of the sample signal auto covariance obtained 
in Chapter 4 have been input to an adaptive filtering program. This program 
uses simulated range rate observations to determine the error and uncertainty 
of the broadcast state vector after processing three new observations at every 
integration step treating the precise ephemeris state vector as errorless for 
this purpose. 

Conclusions and recommendations are made in Chapter 6. 

1 . 2 Brief Description of Data Utilized 

The data utilized in this study was received from several sources, a 
brief description of which is given below. The first data set (hereafter referred 
to as Data Set L) consists ol ptveise and NAC predicted state vectors for 


i 


satellite nos. 13 and 19 for several passes during a period in September, 1974. 
This data was received from Defense Mapping Agency Aerospace Center 
(DMAAC). The data is in the earth-fixed coordinate system of NWL 9D as 
described in (White et al. , 1975]. As mentioned earlier since 1975 the broad- 
case ephemeris system has been upgraded. Along with changes in computational 
procedures, the coordinate system adopted has been changed to the Department 
of Defense World Geodetic System 1972 (WGS 72) which differs by a small amount 
from the NWL 9D system in which the precise ephemeris continues to be main- 
tained. 

The second data set, hereafter referred to as Data Set D, was received 
from DBA Systems, Inc. It consists of Doppler observational data from three 
ground stations in Florida acquired during a period in January, 1976. As JMR 1 
receivers were used, the broadcast ephemeris for the satellites tracked was 
also available in a message form. The precise ephemeris of satellite nos. 19 
and 20 for the related period was obtained from DMA . 

The third data set, hereafter referred to as Data Set S, was received 
from the National Geodetic Survey (NGS). It consists of the precise ephemeris 
for satellite nos. 12 and 19 tracked during October, 1976. The broadcase ephem- 
eris of the satellites for the related period was received separately from the NAG. 

1 . 3 Brief Description of Computer Software Utilized 

This study has required a considerable amount of data processing for 
which software from the following sources was used after due modifications. 

For obtaining station positions with Doppler observational data, the 
Short Arc Geodetic Adjustment Program (SAGA) as received at The Ohio State 
University and described in (Kumar, 1976] was used along with the stand-alone 
program SAMVAP received from Air Force Geophysics Laboratory. For de- 
coding the JMR 1 receiver ephemeris message, a routine obtained from Mr. 

White, DMAAC, was useful. Fbr Kalman filtering procedures, the program 
"KARTHOD" from the University of Texas at Austin was suitably modified. 
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1.4 Coordinate Systems Used 


For formulations related to satellite dynamics, an earth- centered 
inertial coordinate system (ECI) has been used as defined below: 

The X and Z axes are directed, respectively, to the true vernal equi- 
nox and the true North celestial pole at a selected epoch t Q . The V axis forms 
a right-handed system with Z and X. 

The precise and broadcast ephemerides provide the satellite state 
vectors in an earth-fixed system (EF), through there is a small difference in 
the scale and longitude definition between the two. The earth-fixed coordinate 
system is defined with the X axis oriented to the Greenwich Mean Astronomical 
Meridian, and the Z axis passing through the Conventional International Origin 
(CIO), both as defined by the Bureau International de l'lleure (BIH). The Y 
axis forming the right-handed system defines with X the average geodetic 
equator. 

For observations, topocentric systems which are parallel to the above 
but passing through the observer position instead of being geocentric are also 
used. 

Variations from the above where they arise, and actual symbols used, 
have been explained in the text. All vector quantities have an overbar. 
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2. RANK DEFICIENCY PROBLEM IN THE DOPPLER SYSTEM 


2. 1 Introductory Remarks 

One of the practical aims of geodesy is that of the determination of 
positions of points on the earth's surface. The latter aim has dictated the 
need for adopting a frame of reference (or coordinate system) with respect 
to which locations of points could be determined. Having adopted a refer- 
ence frame, it is imperative that for the results to be fully meaningful the 
coordinate system is maintained. 

Unfortunately, one kind of geometric observation cannot provide all the 
necessary information about the coordinate system. For example, range 
observations can give information about the scale but not the origin or orienta- 
tion of a coordinate system . When an adjustment is carried out with such 
observations by the usual method of observation equations [Uotila, 19G7] with 
the mathematical model 

Y a - C(X.) (2.1) 

where G is a vector function relating the u X 1 parameter vector X a with Y a , 
the n x 1 observation vector, the following linearized form is obtained 

V = HX + Y (2.2) 

where 

V is the n x 1 vector of observational residuals 

H is the n x u matrix <3G/dX a , the matrix of partials of the observables 
with respect to the parameters evaluated at the nominal value of the 
parameters Xo 

X is the u x 1 vector of unknown corrections to X 0 
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and 


Y ■■ Yo - Yb is the difference between the computed obseiwations 
Yo Ci (No) and the observations Y b . 

Generally in all geodetic adjustments n > u, and the rank of the H matrix becomes 
less than u if no parameters are constrained. 

So in the normal equations 

NX -U (2 .3) 

win' re 

N ll T Pll, the u x u normal matrix 
P n x n weight matrix of the observations 

U Il'PY 

the rank of N is less than u, and a Cayley inverse cannot be obtained for 
N, which could give the unique solution X -N l U, as the unbiased estimator 
of the parameters. 

This is because the lack o f information about the reference frame in 
the observations leads to a rank defect of the design matrix H and the singu- 
larity of the normal matrix. 

To overcome this situation, additional information about the coordinate 
system lacking in the observations is needed in the form of constraints on 
parameters. This additional information may be introduced in the form of 
(a) inner constraints, (b) weighted constraints, or (c) absolute constraints. 

The nature of the constraints enforced influences the coordinate system in 
which the solution vector is obtained. 

For making decisions in this matter, the rank deficiency encountered 
in the short arc mode of Doppler data reduction used in this study has been 
derived in this chapter and the decisions regarding the constraints have 4 been 
explained. 
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2.2 Rank Defect in the Short Arc Mode 


Primarily, the rank defect in a design matrix depends on 
the type of observations. However, the rank defect may increase due to the 
additional parameters used in the mathematical model other than the station 
positions and/or numerical problems, which may generate dependence between 
the columns of the design matrix. 

Therefore, only the minimum rank defect situation has been analyzed, 
taking into consideration a simplified form of the mathematical model used in 
the SAGA program with the Doppler data reduction [Brown, 1969, 1973], and 
this is termed the fundamental rank deficiency. 

In its simplest form the integrated Doppler shift is modeled as 

D = Af (t 2 - ti) + (p P2 - p Pl )/X (2. 4) 

where 

X is the wave length of the adopted frequency of transmission 

D is the observable Doppler count at one ground station P during 

the motion of a satellite from position 1 to position 2 
t lf t 3 are the epochs corresponding to satellite positions 1 and 2 
Af is the unknown frequency offset given by f e - f 8 , where f g is the 
ground reference frequency and f g is the frequency transmitted 
by the satellite 

Pp 2 *Ppi are the slant ranges corresponding to positions 2 and 1 from P. 

In terms of the computed range difference, the observation equation in 
a functional form is given by 

(Pps - Ppi) + dPi2 (dX P , dXj, (K 2 ) - (Pp 2 - Ppi) + V 12 (2.5) 

where 

(Pp2 - Ppi) c is the computed range difference 
(Pp2 - Ppi)° is the observed range difference 
[D - Af(tc - tx» X 

X [D - A f° (t 2 - tx)] - X (tj. - ti) dAf 
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Pi? Pp? “ Pp 1 

dPi ? is the change in range difference due to changes in coordinates 
of station and satellite positions 

V 12 is the observational residual corresponding to station P and 

satellite positions 1 and 2 
dX P are the three ground station unknowns at P 
dXi arc the three unknowns for satellite position 1 

dXj are tK three unknowns for the satellite position 2 

dAf is the unknown frequency offset correction to the nominal 
value Af° 


For the study of rank deficiency it is necessary to examine the structure 
of the design matrix H, arising from the left member in the observation equation 
of the i*esulting form 

dpiz (cK Pt dXt, dX 2 ) -i X (t? - ti) dAf - 

X [D - Af° (te - ti)] - (Pp? " Ppi) C + V* (2. G) 

and to determine the number of independent columns therein. Fbr doing this, 
the form of II matrix in a range observation will be considered first since the 
formulation for range difference observation follows easily from it. 


2. 2. 1 H Matrix for Range Observation 

The formulation of the II matrix is simplified in the system of topocentric 
right ascension and declination. Lot 


Xr 


Xi 



be the earth-fixed coordinates of ground station P, and 



be the earth-fixed coordinates of the satellite at epoch 1 


If GO anti 6; a re the topocentric right ascension anti declination of the* satellite 
at position 1 , then 
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Xp-Xi 


Yp - Yi = 

Zp - 7j\- 


1 0 x«l rppx cos 61 cos ai" 

0 1 -y. Ra (6 0 i) Ppi cos 6 i sin 0 £i 

-x. y« lJ Lppx sin6 x 

[Krakiwsky et al. , 1967] 


where 


R3 (6gi) the transformation matrix for rotation around the third 

axis given by r cos 6^ sin 0 Gi 0 " 

-sin 6^1 cos 0 s i 0 


x B , y, 


is the Greenwich apparent sidereal time at epoch 1 
are the components of polar motion 


Denoting the polar motion matrix by C T , equat ion (2 . 7) can be differen- 
tiated treating the slant range Pn and the position vectors X P and Xj as vari- 
ables. 

Neglecting second-order terms, C T is orthogonal. Therefore, 


-dXp -dX x 
C dY P - dYi 
-dZp - dZj. 


cos («i - 0 G i) cos 6f 
sin (ai - 0 Gl ) cos 6x dp P i 
sin 61- 


Therefore, 


"dXp - dXi 

[cos (Ox - 0 G1 ) cos 61 ; sin(ai - 0 Gl ) cos6 x ; sin6x] C dYp - dYx 

.dZp — dZi 

This is the familiar observation equation form. 

Denoting 


- dp P i 


dXi = dYi 


equation (2. 8) can be rewritten as 
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T 


T 


(2.9) 


dX P 
(fifi 

where h P1 is a 1 x 3 submatrix given by 

hpi - [cos(a x - 0 Gl ) cosSj + x„ sin6 x ; sin(Otx - 0 Gl ) cos fij - y B sin6x: 
- x B cos(dfi - 0 Gl ) cosfii + y B sin(a x - 6bi) cos 6, + sinfix] 

This equation gives the structure of the design matrix in the case of a slant 
range observation Ppi . 


dp P1 - 


|ji P x -h«j 


2. 2. 2 II Matrix for Range Difference Observations of Doppler System 

To obtain the structure of the design matrix in the case of a range dif- 
ference observation, equation (2.9) is extended to consider corrections to ranges 
dppx and d /0p 8 . Corresponding to (2.9) for dp Pl , 

dXp 

dX 8 

Hence dPx 8 , the variation in range difference due to variation in position of 
station P and satellite positions 1 and 2, is given by 



dPxs 


dPpi - dpp2 


[hpi 

“hpi 

ol 

"dX P " 

dXj 

- r hps 

0 

"bps j 

’dXp" 

dXx 

L 

J 

_dX ? - 

L 


J 

.dX ? . 


c, 


|(hpx - hps) -hpx h 


c 


•] 


'dX P ' 
dXx 
. dX 2 . 


(2.10) 


The variation of the range difference due to variation of the last 
unknown dAf is now considered. This is readily seen from equation (2.6) as 
X (t?-tx). Thus, the structure of the 11 matrix for a Doppler observation for 
one ground station, two satellite position unknowns, and a frequency offset 
unknown will be 
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(V) 

- sin(0(i - 6^) cos 61 


(i) (ii) 

cos(Oti - 0 s i) cos 61 sin(OCi - 0 8l ) cos 61 

- cos(Of 2 - 0(2) cos 62 - sin(a 2 - 0^) co 

(iii) (iv) 

sin 6i - sin 6 2 - cos(Oi - 0 8i ) cos 61 

(vi) (vii) (v iii) 

- sin cos(a 2 - 0 82 ) cos 6 2 sin(a 2 - 0 82 ) cos 6 2 

(be) (x) ] 

sin 6 2 X (tc - ti) 


( 2 . 11 ) 


treating the polar motion components as known and leaving them out of con- 
sideration for rank deficiency analysis. 

Now the addition of a second range difference observation for satellite 
positions 2 and 3 will imply three more unknowns. Denoting oo t = a { - 0 a i, the 
structure of H matrix for the 13 unknowns with two observations will be 

(i) (ii) (iii) (iv) 


cos CO! cos6! 

sinooi cos ^ sin 61 - 

sin6 2 

- COSOOl COS 61 

- COS U) 2 COS 62 

- sin oo 2 cos 62 




cos tOfe cos 6 2 

sin ufe cos 6 2 sin 62 - 

sino 3 

0 


- cos C0 3 cos 63 

- sinoo 3 cos 63 




(V) 

(vi) (vii) 

(viii) 

(ix) 

- sinco! cos 61 

- sin cos co 2 cos 6 2 

sin 002 cos 6 2 

sin 63 

0 

0 -cos 002 cos 6 a 

-sinoo 2 cos 6 2 

- sin 6 2 

(X) 

(xi) (xii) 

(xiii) 



0 

0 0 

*<tc - tl) 


(2. 12) 

COS Uh COS 63 

sinco 3 cos 63 sin 63 

A(to - ts) 



2.2.3 H Matrix in SAGA 

It can readily be seen that the above procedure gives three more unknowns far 
every additional observation, and an overdetermined system required for adjust- 
ment cannot be obtained. This situation is remedied either by increasing the 
number of coobserving stations or restricting the number of satellite unknowns 
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1 


per pass. In SAGA, the satellite unknowns are restricted to six per pass 
(three for position and three for velocity) by assuming the force model to be 
known. In this investigation the same procedure has been followed. 

Adopting the following compact notation 

CC = cos to. cos 6. 

i i i 

SC = sin to. cos 6. 

i 11 

Sj = sin 6. 

the problem is extended to one ground station and four consecutive satellite 
positions. Following equation (2.12), the structure of H matrix for three obser- 
vations and 16 unknowns assuming one frequency offset per pass, will appear as 


(i) 

(ii) 


(iii) 

(iv) 

(v) 

(vi) 

(vii ) 

(viii) 

(ix) 

cc x -cc s 

sc 5 -sc 3 

Si 


-CC X 

-sc x 

-Si 

cc s 

sc 2 

S8 

CC 2 -CC3 

sc 2 -sc 3 

Ss 

-sb 

0 

0 

0 

-CQ 

"" SCg 

“S3 

CC3- cc 4 

sc 3 -sc 4 

Sb 

-S 4 

0 

0 

0 

0 

0 

0 


— dXp 

— 

—i 

( 

- dxi — 

— ) 

( 

■ dX 2 “ 

— 



The unknowns for satellite positions are now reduced to six, viz., dX 0 , 
dX Q , the corrections to assumed values of position :md velocity components 
at an adopted epoch to (preferably taken as the epoch t 0 at midarc of the pass) . 
This approach gives the following linear transformation [Brown, 1969, p. 20] 
for satellite position i at epoch tj: 

rcero- 

dX, R 1 ( 2 . 11 ) 

:n 33 33 vrr 

L dXo J 


U 
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where [X 0 T X 0 T 1 [Xo Yo Zo Xo Yo Zo]» and O* is the matrizant in the 
inertial system given bv 

d (X|, Yi, Z|, T)i 

O 1 - . , 

d (Xo , Yo , Zo , Xo , Yo , Zo ) 

where 

Xi, Y|,Z i are the geocentric inertial coordinates of the satellite at 
T, with Ti = tt - to 

Xo , Yo , Zo , Xo , Yo , Zo are the assumed initial conditions at T 0 
The inertial coordinate system is defined as the coordinate system which 
is coincident with the earth-fixed system at T 0. 


O is obtained in SAGA with the help of the orbital integrator developed by 
Hartwell fl9(>8|. it employs a power series solution to the equations of motion 
in the inertial system, in a recursive algorithm, which can be represented by 
the following for an arbitrary Tt: 


Xi xr 




a? 


•V 

Yi Y, 


bo 

bi 

br 


u, 

-7.i Z,- 

1 

-Co 


Co 

• . • 

Ca- 


LT‘ , qr* U 

i 


(2. 15) 


where the coefficients a,, b„ e Jt i 0, 1, .... q are functions of the six 
initial conditions at T 0, (X,,, Yo. Zo, Xo, Yo, Zo) and the earth gravitation:) 1 
coefficients, q is the index of the power series at which the series is truncated 
when a prespecified tolerance (0.001 m) is reached for T„ the maximum value 
of T to be exercised. 

From (2. IS), 0 can be seen as 



— i 

p 


i — 

c 

n 1 

n..i 

< 

a 



-0:« 

O.e 



where is, in turn, a polynomial 1 l, 9; n 1, (’». 


If. 



(o^ Gfi 


Oln 




1 

r 


T 


4 


( 2 . 10 ) 


- i 

wit luOfci as the dominant term. The correspondence between the a coefficients 
of (2. U'O and (a,b t cl coefficients of (2. 15) is easily seen. For example, taking 
On. (t>o)n 


K 


t 


" cosv t i sin V T t 0 ~ 

- sin V Ti eos0 T { 0 

0 0 I - 


(2. 17) 


is the matrix which can transform the matri/ant O 1 to an earth -fixed system 
at ti. V is the earth rotation rate. 

For passes up to 20 minutes, T< ; 10 minutes, cose T j 1.0, and 

sin V Ti * 0, so R 1 is taken as an identity matrix for the purpose of rank 
deficiency analysis. Thus with the help of expiation (2. 1 I), expiation (2. 12) is 
used to obtain the corresponding portion of the H matrix for satellite position i 
in terms of the new satellite unknowns through the linear transformation 


r i 


MX/ 

h ; 

[<ix t 

h. n‘ 

9 

1 :i 

a l 

l 3 o 

.. il\ 0 . 


a i 


For the situation in (2. 12), the unknowns will reduce from Id to 10, 
and the new 11 matrix will have the structure 


Id 



(i) 


(id 


(iii) 


(iv) 


cc ] -cc 2 ; sCj-sCg Isj-Sg :-(cc 1 )frf 11 )-(sc 1 )(0 1 21 )-(s l )( d 3 \) . 

| i ;+(cc 2 )( n 2 11 )+(sc 2 )( o 2 2 i) (s 2 )( o ? 31 ) ! 

1 , 11 I 

i i 

CC 2 -CC 3 | SC 2 -SC 3 j S 2 -S 3 : r (CC 2 )(0? 1 j -(SC 2 )/ rf 21 )-(S 2 )f C?31) j 

• ! J+(cc 3 )(j? ll ) +(sc 3 )in? 2] ) +(s 3 )( 0 3 3 i) : 

1 ! 1 

1 : 

cc 3 -cc 4 1 sc 3 -sc 4 ;s 3 -s 4 ; -(cc 3 )(rfn)-(sp s x ^ 2 i)'< s 3 )< £ ^ 3 i ) i 

: : :+(cc 4 H.n 4 u )+(sc 4 )( 0 4 21 )+(s 4 KQ 4 31 ) 

(v) (vi) 

! -(cc 1 )( 0 1 12 )-(sc 1 )(d 22 )-(s 1 )(«! 32 ) i-(cc 1 Hn 1 13 )-(sc 1 )(oi 23 )-(s 1 K«fi! 33 ) 

■ +(cc 2 )(fi 2 12 >f(sc 2 Hr? 22 >f(s 2 )(rf 32 )|+(cc 2 )/o 2 13 )+(sc 2 )(d 23 )+(S 2 )(rf 33 ) ! 


' -< cc 2 )(Q 2 i 2 )-< s ‘ c 2 ) ( n 2 22 )-( s 2 Hd 32 ); -(CC 2 )( 0 2 13 )-(SC 2 )(rf 23 )-(S 2 H n 2 33 ) • 

:+ (cc 3 )(n 12 )+(sc 3 )f n 3 22 )+(s 3 h ft 3 32 ). +(cc 3 )( 0 3 13 )+(sc 3 )id } 23 ) + (s 3 )(o 3 2B ) j 


; ) 

; -(CC 3 ) (d 12 )-(SC 3 )((i 3 22)-(S 3 )<c? 32 ) : ~< CC 3 )(d ls ) "(SC 3 )( jf 23 )-(S 3 ) ( C?33> 


A ° I 

: + (CC 4 )( 0 12) + (SC4) C2 22 MS 4 )(0 3 2 \ + (OC4)(f? 4 i3) + (SC 4 )(^ 2 ^j , d 33 ) • 

(vii) (viii) 

! -(CCpm^l-fSCjK nl 24 )-(S 1 )(rt 34 , :-(CC,)(» 5 HSC 1 Hrt 25 )-<S 1 M n l 35 > 

i + (CC 2 )(f^| 4 >* (SC 2 ) rf 2 4 |4 ‘< S 2 l,n 34* ,CC 2 I < n2 15t*(SC 2 )( ! * 25 )+(S 2 )(f* 35 ) 


•(CC 2 H£# 14l -^c 2 )(rf 24 )- 6 2 )(02 34 ) !-(CC 2 xn 2 15 ) - (SC )( t# ,, 

' J 2 25'~< & 2 ,< 


rf. 


+ ( CC 3 


, 3> ' ° ° V<V< ^34* *<«V< a3 ^ ){ 

• ! 1 
;' ,cc 3 ,( n* m’' (SL 3 i( n 3 24'-<s 3 »n 3 34 )i-(cc 3 K c^ 5 >-<sc 3 Hn 3 25 l - ( s 3 i<o 3 15 )! 

i+icc 4 k ft 4 i 4 i*(sc 4 i(n 4 24 x(s 4 K ^ 34 ) ;mcc 4 x n 4 | 5 n ,sc 4 ), n 4 25 n ,s 4 )( n* 
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-(CCjWfi 1 jgWSCjHfl^WSjX Q^g) 
+ (CC 2 )l « 2 16 )4 (SC 2 )({f *■ (S 2 )( 0 2 3 g) 


^2"V 


I -<«V ( °V<sc 2 x n 2 26 )-.s 2 )( o 2 3( .i 

; MCCjM n3 16 t*(SC 3 K n 3 2 6 )t ( s 3 l( rf I 


, X W 


! -< cc 3 Hn i 1( .H S c 3 „n 3 2 6 )-(s 3 i(0 3 36 , 


X (t 4 -t ) 


Using index u for the satellite position and v for the column of the matrizant 


C2, and denoting 

CC + SC + S by t 

u lv u 2v u 3v 

and (t - t ) by T 

' u* l u u+1 , u 

the above matrix takes the following more compact form 


m tr 

u v 


<i> 

(>») 

(iii) 


(iv) 

(v) 

CC-CC 1 SC -SC 1 
1 2, 1 2, 

1 ' 

s r s 2 

-ml 

1 


i -m, (i n 4 ni fl? 
,12 2 2 

! 

CC -CC„ 
2 3 

, 1 
isc -sc ' 
, 2 J l 

S 2’ S 3 

: - m 2 

^l 4m 3 ^1 

! 

: -n, 2 rf 2 .m 3 C? 2 

j 1 

_c r 3 -cc 4 , S c 3 - S c 4 ; 

S 3" S 4 

i 

i 

-m 

! 3 

3 n 

Oj+m4“ j 

'.- n, 34 +m 4 n4 2 
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(Vi) 


(vii) 


( viii ) 


(ix) 


:- m . n V m 2 rf 3 : • m i rf 5 +m 2 c# 5 


I 


-m 




i 2 3 3 3 


|-m 3 rf 3 *m 4 n 4 3 






I 


2**4" ul 3 4 ‘4 ! - m 2 ^5 rm 3* # 5 J -« V 6 +II1 3 * ' 6 


-m l d 6+ m 2 Q 6 


rf c nn,rfc i ft 3 . 


-m 


C? .+m A d 


.4 i 


3*'4 “V*4 :' m 3"5 ‘T5 ; ‘“3*'6'“‘4”6 


rf +m i f^„ i -m„Q 3 +m , 


) (_ 


dX 


(x) 
X r 


21 


X T 


32 


X r 


43 


dAf 


2.2.4 Rank Defect Analysis 


( 2 . 20 ) 


As further observations are taken in the same pass from the same 
station P, the number of unknowns do not increase, and only more rows will 
be added to the above H matrix (2.20) which can now be analyzed for rank 
deficiency. 

Examining the first row, it can be seen that 
element h 14 = -mi fix + ma Of 

= -(cos u>i cos &i) On - (sinc^ cos 6i) fia - (sinSofla 
+ (cosco 2 cos^s) ftu + (sinOJa cosS^fia + (sin6a)f23! 

- - hn ~ hiaQa - hi3 fill (2.21) 


PAGE 13 

Ob r\TT a * 


19 


as On — flfx* Ojj — flgx. ^31 - flax, flix ~ fl?i 3 s they are monotone 
decreasing series with coefficients which alternate in sign [Hartwell and 
Lewis, 1967] and have the same first dominant term (a 0 ) u independent of T 
as seen from (2.16). Similarly, fl a - fl a and flax - fllx. as they have the 
same first dominant terms (<*o)a and (O' 0 ) 31 , respectively . 

Therefore, hx 4 is a linear combination of hxx, hx 2 , and hxa. Similarly, 
it can be shown that his, h 16 , hi?, hie. and hxsare linear combinations of hxx, hi 2 , 
and hx3. And this holds for everyrow of H (e. g., ha*— -h a O?i - hasflfx - l^flli). 

Also, from row 1 to row 2, the scalars in the linear combination will be approx- 
imately equal in a short arc, e.g„ flnsa flfi, as explained earlier. The same argu- 
ment holds from one row to the next, and columns 4, 5, 6, 7, 8, and 9 become 
dependent columns leaving only four independent columns . 

If more stations observe the same pass, there will be three more station 
unknowns per station, but the above arguments about rank deficiency will hold. 

It can therefore be concluded that the fundamental rank deficiency in a Doppler 
system, short arc mode, is 6. 

It is assumed that the scale information is obtainable in the system 
from A , the wave length of transmission, as seen from the tenth column of 
the H matrix above. Having determined the fundamental rank deficiency in 
the system, an effort will now be made to determine what quantities are esti- 
mable in the above situation. But it is obvious that a rank deficient matrix 
like H above will lead to a singular normal matrix N. 

2.3 Estimable Quantities in Doppler Obsei’vations 

It is known that if a normal matrix N is singular and a solution is obtained 
with a pseudo- inverse N + , the solution vector X. - - N U is not estimable since 
E(X a ) / X, the parameter vector. 

As derived by Rao [1973], any arbitrary matrix G can make GX, estimable, 
in a linear system, if the condition 

G[I-( N) + ( N)] - 0 (2.22) 
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is satisfied, where I is the identity matrix. 

So if a change can be made in the parametrization of the mathematical 
model obtaining an H matrix which is not rank deficient, ( N ) = ( N )*, and 
the new parameters will be estimable. This can be done in the Doppler system 
considering the form of the H matrix given by (2 . 20) . 

With the modified parameters 



instead of X p and X 0 , the restructured form of H will be as given below: 


CC f CC 2 +m l fl r^^l ! SCj-SCg+njj O 1 -I HgCfg i S fV nij rfg- 


CC 2' CC 3 + “2 °V ^1 | SC 2 _SC 3 + m 2° 2 2 -m 3° 3 2 . VV 


.2 _3 


. CC 3" CC 4 +m 3° 3 l" m 4° 4 l I SC 3" SC 4 +m 3 Q3 2 -n M^2 ' W m 3 °V ^ 


dAx^ 
po 

( IV) (v) (vi) (vi 

I _ 1 2 1 2 ' 1 2 ' 

! 1^4 4 : " m 1^5 +m 2^5 | _m 1^6 +n ^2^6 | * T 21 


' _ 2 3 i 2 3 ! 2 3 I 

4 “b 0 4 1 " m 2^5 +m 3 0 5 i ” m 2°6 +m 3 n 6 ' ^ T 32 


‘ 3 . 4 


! m3 0 4 4 'pi 4 0 4 | * m 30 5 +m 4^ 5 , ” n b^6 +m 4^6^ 


43 I 
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In this matrix all the columns are independent, and the solution will be a vector 
of estimable quantities. Thus, it is concluded that in a Doppler system, the 
vector of coordinate differences between the observing station and the mid arc 
state vector of the pass, the velocity components of this state vector, and the 
frequency offset are estimable. 


If more than one station is coobserving the same pass, the linear 
relation 


AXpp 


-Xp- 


*X Q - 


-X P -Xo- 


-Xp-Xo-j 

Yp 


Yq 

= 

Y P - Y 0 

- 

O 

>* 

t 

£ 

-Zp- 


-Zq- 


- Zp - Zq- 


-Z q ~ ZqJ 



can be used to conclude that interstation coordinate differences are estimable 
if the stations coobserve the same pass of the satellite in a short arc mode. 
The coordinate differences are independent of origin. The scale information 
comes from the wave length of transmission, and the orientation information 
comes from the force components enforced in the satellite dynamics. 



3 . EPHE ME RIDES O F THE NNSS AND 
THEIR ACCURACY ESTIMATES 

3.1 Introductory Remarks 

Having analyzed the rank deficiency problem in the previous chapter, 
the next step in the goal to achieve improved positioning is to assess the accu- 
racy of the ephemerides of the NNSS. As mentioned in Section 1.1, the 
ephemeris of satellites of the NNSS are available in two forms, precise 
ephemeris computed after the fact by DMA and broadcast ephemeris which 
is obtainable on a real time basis from satellite transmissions . 

These values can be treated as direct observations on the position and 
velocity of the satellites for applying conventional sampling techniques to 
obtain estimates of uncertainties . Estimates of precision or more correctly 
the prediction errors in the broadcast ephemeris can be found by comparing 
the two values of state vectors for common time points in the overlaps of 
successive orbit fits. 

Since the precise ephemeris is known to be more accurate than the 
broadcast ephemeris, in the pursuit of improving the broadcast ephemeris, 
estimates of its accuracy have been computed by comparing it with the precise 
ephemeris. After a review about the satellite system and ephemerides, the 
results of this study will be presented. 

3.2 The Navy Navigation Satellites 

There are at present (May, 1977) five operational NNSS satellites in 
orbit. Their identification numbers are summarized in Table 3-1. Some 
typical orbital elements of the satellites as per the latest data set (Data Set S) 
are given in Table 3-2. The orbital elements tabulated are the mean motion (n). 
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Table 3-1 

NNSS Satellite Identification Numbers 


Launch IJute 

Apr 14,1867 

1 

h* 

% 

Sep £5,1007 

Mar 2, 1006 

Aug 27, 1070 

Oct 29, 1073 

Preciae 
Kphemet'ia 
Utmlf float ion 
Number 

58 

50 

60 

♦ 

6 ! 

68 

77 

Operational 
Ephemerl* 
Ident if i cation 
Number 

30120 

30130 

3014 0 

30180 

30100 

30 2 0 0 

Number 

1 






Used in 
lbie Stiuty 

12 

13 

14 

18 

19 

20 


* SoiellUe Number 18 haa <iinee been declin ed a# non operational. But tome data far thta 
rt itullite was available in thia niudy in Ihtta Bet P, 


Table 3-2 

Values for Selected Operational Ephemeris Parameters 
for Selected Days in 1976 


Satellite 

Hay 

n(deg/mini 

e 

n(km) 

CO# i 

12 

329 

3 . 3616808 

. 002166 

7440.73 

>0.004096 

13 

316 

3 . 3663808 

. 001076 

1 1 

0.006710 

14 

316 

3.3720148 

. 004000 

7463.47 

0.013237 

19 

320 

3.3657246 

.017600 

7464.09 

-0. 002161 

20 

330 

3.4 100666 

. 

.016762 

7306.34 

* 0 . 002291 
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eccentricity (e), semi-major axis (a), and the cosine of inclination (cos i). 

All NNSS satellites follow a near circular polar orbit. 

3.3 Precise Ephemeris 

Precise ephemerides for one or more Navy Navigation Satellites are 
conputed on alternate days based on 48 hours of observation made at over 20 
stations distributed around the world [Anderle, 1976]. 

The equations of motion of the satellite and the variational equations for 
the forces are numerically integrated by a tenth-order Cowell process with UTC 
time as the argument of integration. The force equation includes terms for the 
gravitational field of the earth, moon, and sun, the lunar and solid earth tide 
effects, atmospheric drag, and solar radiation pressure. The gravitational 
field of the earth is given in a spherical harmonic expansion containing about 
400 terms, and the earth tides are based on a Love's number of 0.26. The 
gravity field revised in January, 1973 (NWL-10E) is in current usage. 

The precise ephemeris is believed to have periodic errors of about 2m 
in each coordinate due to uncertainties in the earth's gravitational field and 
effects of variations in atmospheric density on the computed satellite positions. 

It is maintained in the NWL-9D coordinate system which is believed to be related 
to the NWL-10F system (consistent with WGS 72) as follows: 

longitude Xi 0 r - X 8D + 0'.'260, X is positive East 

geocentric latitude 'I'lor ^ 'I'bd 

radius £ 10 r ~ £»o - 5.27 m 

[Anderle, 1976] 

Based on the above approximate relations Vincenty [1976] has derived the 
transformation parameters given in section 3.5.2, which have been used in 
this investigation. 

The precise ephemeris is made available at one-minute intervals in the 
geocentric earth- fixed system in the form of position and velocity components. 
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3.4 Broadcast Ephemeris 

The broadcast ephemeris is computed as already explained in Chapter 1. 
Each satellite has a memory which can hold 16 hours of orbit prediction data. 

This predicted ephemeris is injected into the satellite memory about every 
twelve hours. 

The broadcast ephemeris is received at an observing station in the form 
of coded parameters from which earth-fixed satellite positions can be calculated 
[Moffett, 1973]. These parameters are divided into 14 fixed orbit parameters 
whose values change only twice a day and four sets of variable orbit parameters 
whose values change every two minutes. These are listed in Tables 3-3 and 
3-4 as taken from [Wells, 1974] and [Moffett, 1973]. 

The decoding of the parameters and the computation of the positions of 
the satellites in the earth-fixed coordinate system, at two-minute intervals, is 
done according to procedures described in [Moffett, 1973]. For velocities of 
the satellite, at two-minute intervals, time derivatives of the variable paramet- 
ers are also required, and these have been obtained by a polynoroinal fit to a 
maximum of 16 consecutive values. 

Since December, 1975, the broadcast ephemeris system has been up- 
graded [Black, 1976], some main features of which are given below: 

(i) The previous A PL 4.5 geopotential model has been replaced by the 
WGS 72 model. 

(ii) The value of GM has been changed from 398 601. 5 ± 0. 6 km 3 /sec 3 
to 398 600.8 ±0.4 km 3 /sec a . 

(iii) Station coordinates of tracking stations have been changed by small 
amounts to bring greater internal consistency. 

(iv) Implementing the main sun-moon-induced body tide corrections. 

The (single pass) error budget is given in Table 3-5 both "before" and 
"after" the introduction of WGS 72. It is clear that the error budget will continue 
to be dominated by uncertainty in the geopotential model and incorrectly modeled 
surface forces. 
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Table 3-3 Fixed Orbit Parameters 


Symbol 

Definition 

Broadcast 
Units 
(Current 
Resolution ) * 

tp 

Time of first satellite perigee after last 
satellite injection 

10' u min UT 

n 

Mean motion (only fractional part is 
broadcast ) 

1 

.7 

10 deg/min 

w(tp) 

Argument of perigee at tp 

10" u deg 

l»l 

Absolute value of precession rate of perigee 

_7 

10 deg/min 

e 

Eccentricity of orbit ellipse 

10” 6 

a 

Mean semi -major axis of orbit ellipse 

10 metres 

n(tp) 

Right ascension of ascending node at tp 

10 _U deg 

e 

n 

Precession rate of ascending node 

10” 7 deg/min 

cos i 

1 

Cosine of inclination 

| 

ID" 6 

GAST(tp) 

Greenwich apparent sidereal time at tp 

\0' U deg 

- 

Satellite identification number 

- 

- 

Day number and time of last satellite data 
injection 

2 min UT 

sin i 

Sine of inclination 

io- 6 

m 

Fractional satellite frequency offset 

(f * f s> /f o 
o so 

parts in 10^ 


* For each of these parameters there la a trailing aero digit which is not 
currently used, and which could be used to increase the resolution (the 
frequency offset has four trailing zeroes). 

Note: The above table is based on Wells [1974]. 
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Table 3-4 


Variable Orbit Parameters 


HSSim 

1 

Definition 

Broadcast 

Units 

t 

Time in even minutes of UT, modulus one half 
hour 

2 min UT 

AE(t) 

Correction to eccentric anomaly at time t 

10’ U deg 

Aa(t) 

Correction to semi-major axis at time t 

10 metres 

T>(t ) 

Out of plane orbit component at time t* 

10 metres 


* »l(t) values are available only at four minute intervals (for times 
which when expressed in minutes UT are divisible by A). 

AE(t ) and Aa(t) values are available at two minute intervals (for 
even minutes UT). 

Note : The above table is based on Wells [1974]. 
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Table 3-5 NNSS Error Budget * (Single Pass) 


Meters 

1 

Uncorrected propagation effects (3rd order 
ionospheric and neglected tropospheric 
eff'cts ) 


1-5 


2 

Instrumentation (oscillator phase jitter) 


1-6 


3 

Uncertainty in the geopotential model 

1 5 - 2 0 ( A PL 4 . 5) 


5- 10(WGS - 72) 

4 

Incorrectly modeled surface forces 
(secular error growth due to incorrect 
period, drag and radiation pressure) 


10-25 


5 

Unmodeled UTI-UTC effects and in- 
correct coordinates of the pole 


* 

X 


6 

Ephemeris rounding error (last digit of 
ephemeris is rounded) 


5 



Overall Uncertainty 

19-3 3m(APL4.5) 


12-2 8m(WGS-7 2 ) 


Adapted from Staff of Applied Physics Laboratory (1975] 



3.5 Accuracy Estimates 


Accuracy estimates of the precise ephemeris have been mentioned in 
Section 3.3 above, as available in the literature. Estimates of precision of 
the broadcast ephemeris can be made by carrying out a comparison between 
the broadcast ephemerides based on two successive orbit fits, in the overlap 
period, at common epochs. 

Estimates of accuracy of the broadcast ephemeris are possible to be 
found by comparison between the precise and broadcast ephemeris at common 
two-minute epochs if the precise ephemeris is considered errorless for this 
purpose. 

3.5.1 Precision of Broadcast Ephemeris 

Ephemerides of satellites for common epochs for the overlap portion of 
two successive orbital fits were available in two forms: 

(i) from Data Set D, from passes tracked during injection 

(ii) from Data Set S, where complete injection information was available 
on tape, in coded form . 

The comparisons between the state vectors of the satellites have been 
carried out in two forms: 

(i) comparison between the position and velocity components in the 
Cartesian system (X, Y, Z, X, Y, Z) 

(ii) comparison between the position and velocity components in 

# # I 

the polar coordinate system (X, <p, r, X, <p, r) 

Since the satellite orbits are polar and near circular, the comparisons 
at (ii) yield results which would correspond very closely to out -of -plane, in-track 
and radial differences. 
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For the above computations, values for X, Y, Z and X, Y, Z were 
obtained from the ephemeris message as indicated in Section 3.4. From 
these, the corresponding X , <p , r, X , <p, r values are readily derived from 
the relations: 

r = (X 2 + Y 2 + Z 2 )^ 

tan,f ’ * (x 2 77 ? <31) 

tanX = £ 

and 


-x- 


r coscp sinX - r sinv> cos X cos<p cosX‘ 


-x- 

Y 

= 

rcos<pcosX -r simp sinX cos <p sinX 


<P 

_z. 


0 r cos <p sinV 


• 

L r J 


which yields 


X 


-sinX 
r cos <p 

cosX 
r cos <p 

0 


" . - 
X 

<p 

= 

- sin<p cos X 
r 

- si nip sinX 
r 

cos<p 

r 


Y 

T 


cos<pcosX 

cos<p sinX 

sin<p 


Z 


t 


(3.3) 


when the transformation matrix is nonsingular. 

Let M be the quantity whose precision estimate is being obtained. 
AM t = M» t - M lt , where M lt is the value of quantity M from broadcast 
ephemeris at epoch t t as per later orbit fit; M, t is the value at the same 
epoch ti as per earlier orbit fit. 


AM 


n 

mean value of M 


(3.4) 


where n is the number of data points. 
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RMS 


Root Mean Square value 


(3.5) 





SD 


£(A Mj - A M) a 


j=i 


n - 1 


Standard Deviation 


(3.6) 


i « « * * • 

The above values for X, Y, Z, X, Y, Z, X, <p, r, and X, o, r are given in 
Tables 3-6 to 3-9. 

Based on the preceding results, the following observations can be made 
giving more weight to the values obtained in Data Set S, where the number of 
data points are over 300: 


(i) Taking conservative estimates (based on maximum values), the 
internal consistency of broadcast ephemeris for positions can be taken as about 
10 m, 11 m, and 15 m in the X, Y, and Z directions, respectively. But since 
X, Y, Z coordinates are correlated, it is more appropriate to consider the 
inconsistencies in the in-track, out-of-plane, and radial directions. The 
estimates for these are 19 m, 14 m, and 4 m in the in-track, cross-track, and 
radial directions, respectively. 

Fluctuations to the above extent show no indication of a bias since the 
RMS values are very close to the values of standard deviation. 

The above values have been arrived at after excluding from consideration 
the values for satellite 20 in Data Set S. The unusually high values for satellite 
20 are due to the fact that, as ascertained through private communication, NAG 
has been experiencing periodic fluctuations in the dynamics of this particular 
satellite due possibly to its low perigee height. It does not, therefore, repre- 
sent the general behavior of the Navy Navigation Satellites. 

(ii) The corresponding conservative estimates, for internal consistency 
in velocity components may be taken as 0.2 m/sec each in the X, V, and Z 
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Table 3-6 


Precision Estimates of Broadcast Ephemeris for Position Components in Cartesian System 


Dau 

Satellite 

No of 

Data 


X 

Y 

z 1 

Sat 

No 

Poi nts 

AX 

RMS 

SD 

A Y 

RMS 

SD 

4 Z 

RMS 

SD 




< m ) 

( m 1 

( m ) 

On) 

( m I 

<m) 

(m) 

(m) 

<m) 


12 

20 

0.8 

6. 6 

6.8 

-8.0 

11.4 

8.4 

-12.9 

16.4 

10. 6 


13 

15 

1.5 

7.4 

7. 5 

0. 9 

1 • 5 

1. 2 

0 . 0 

3.2 

3 j 3 

D 

14 

10 

4.6 

NS 

o 

© 

9.3 

5. 9 

7.7 

6. 2 

11.5 

14.0 

8. 3 


19 

15 

3.7 

7.5 

6.7 

0.3 

1.7 

1. 7 

- 1.4 

3-4 

3.2 


* 20 

15 

3. 0 

7. 0 

6. 6 

-9.5 

10.2 

3.7 

-21.0 

22.4 

8. 1 


12 

373 

-0.4 

7. 5 

7. 5 

-1.2 

7.7 

7. 6 

- 0.2 

8. 2 

8. 2 


13 

378 

0. 3 

6. 5 

6. 5 

0. 2 

7.7 

7. 7 

0.8 

10. 1 

10. 1 

S 

14 

354 

0.2 

5 . 2 

5. 2 

a s 

5. 6 

5.6 

0. 6 

5. 5 

5. 5 


19 

421 

O. 4 

9. 3 

9. 3 

a 4 

10. 2 

10.2 

- 0. 4 

14 0 

14.1 


• 20 

416 

%• 1 

18 . 0 

18. 0 

2. 1 

20. 8 

20. 7 

0. 8 

2 8*8 

28.8 


• cxctpUontl vide remarks in section 3.5.1 



Table 3-7 

Precision Estimates of Broadcast Ephemeris for Velocity Components in Cartesian System 


bat* 
S# t 

Saulhu 

No. 

No of 

LftU 

Point# 

x, 

( «n/#j 

V. 

(tn/m 

z 

(ro/#l 

AX 

HMS 

Si) 

IT 

RMS 

SD 


RMS 

SD 


12 

20 

. 028 

. 033 

.0 16 

-.010 

. 036 

. 036 

- 0 18 

. 08 3 

. 083 


13 

IS 

. 0 5 4 

. 056 

. 0 14 

. 004 

. 038 

. 039 

- 018 

. 080 

. 08 1 

D 

14 

10 

. 024 

. 040 

. 034 

- 021 

. 04 0 

. 036 

024 

. 073 

. 073 


13 

IS 

.029 

. 039 

- 027 

002 

. 046 

. 04 6 

- . 030 

. 100 

098 


20 

15 

. 040 

. 045 

. 021. 

. 002 

. 056 

. 058 

- - 010 

. 069 

. 07 1 


12 

373 

00 1 

. 0 95 

. 096 

-.011 

.133 

.133 

- . 003 

. 154 

. 154 


13 

379 

. 00 4 

. 136 

.136 

-.0 10 

. 160 

. 160 

-. 009 

. 188 

. 18 8 

S 

14 

354 

. 003 

. 151 

. 151 

-.012 

. 157 

. 157 

013 

. 202 

. 20 2 


19 

421 

. 000 

. 107 

. 107 


. 152 

. 152 

-.018 

. 187 

. 186 


20 

4 16 

-.00 7 

. 14 7 

. 147 

-.019 

. 170 

. 169 

.012 

. 223 

. 223 
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Table 3-8 


Precision Estimates of Broadcast Ephemeris for Position Components in Polar System 


CO 

cn 


HtU 

Set 

Satellite 

No. 

No. o( 
DaU 
Points 

A 

¥> 

r 

2TA 
( ” ) 

R MS 
( " > 

S D 
( ” ) 

w 
< " > 

RMS 
( M ) 

S D 
( " ) 

SF 
( m ) 

RMS 
( m ) 

S D 
( m ) 


12 

20 

. 026 

. 224 

. 229 

-.429 

. 543 

. 341 

- . 6 

4.4 

4.5 


13 

15 

. 053 

. 232 

. 234 

. 007 

. 095 

. 098 

- . 2 

. 4 

. 3 

D 

14 

10 

. 189 

. 3 20 

. 272 

. 385 

. 450 

. 245 

- . 4 

3. 0 

3.2 


lit 

15 

.13 8 

. 272 

. 243 

- . 048 

. 113 

. 106 

- . 0 

.4 

• 4 


20 

15 

. 049 

.192 

. 192 

- .650 

.692 

. 246 

-2.5 

4. 5 

3.9 


12 

373 

-.021 

. 384 

. 384 

- . 020 

. 306 

. 306 

, 1 

2.9 

2 9 


13 

378 

. 044 

i 1 

1 . 320 

.317 

. 020 

. 356 

. 355 

. 0 

3. 2 

3. 2 

S 

14 

354 

.056 

I 

| . 337 

1 

. 333 

.014 

. 19 3 

. 193 

. 0 

2. 8 

2. 8 


19 

421 

.070 

! .319 

.312 

-. 005 

. 502 

. 503 

. 2 

3.6 

3.6 


• 20 

416 

. 000 

.324 

. 325 . 

. 042 

1.068 

1. 088 

. 2 

3.3 

3. 3 


* exceptional vltle remarfcs In aecllon 3.5. 1 




RMS 


(tn/m) 


~ . 037 
- .012 
- . 006 
- . 050 
.000 


169 

126 

095 

227 

076 


. 029 
. 009 
025 
* 029 


. 269 
. 342 
. 361 
• 385 


SD 

EH 

KBS 

SD 

. 169 

- . 022 

. 069 

. 067 

. 132 

- . 021 

. 077 

077 

. 100 

m 

■ 

. 072 

. 227 


. 076 

. 073 

. 078 

028 

. 086 

. 084 

. 286 

001 

.194 

. 194 

. 342 

. 001 

.241 

. 241 

. 362 

. 000 

. 257 

. 257 

* 384 

003 

.218 

.218 


' 016 


424 


002 


. 270 






























directions. The corresponding in-track, cross-track and radial components 
are 0.2, 0.4, and 0.3 m/sec, respectively. 

(iii) At this stage it is also clarified that though estimates of precision 
have been arrived at by comparing ophemerides in the overlap period between 
successive injections, these values more appropriately represent the error of 
prediction over the inter-injection period. This point will be brought up again in 
a subsequent section while discussing the overall accuracy estimate of the 
broadcast ephemeris. 

3.5.2 Accuracy Estimate of Broadcast Ephemeris 

Estimates of accuracy of broadcast ephemeris have been obtained in two 
ways in the cases of Data Sets L, D and S by comparison of precise and broad- 
cast ephemeris of the satellite for which a precise ephemeris was being main- 
tained for (he related period. 

In the case of Data Set L, both the precise and predicted ephemerides 
(broadcast ephemeris prior to injection) were available in the form of earth- 
fixed position and velocity components (X, Y, Z, X, Y, Z) in the same 
coordinate system (NWL 9D); and as explained in (White, 1975] the predicted 
ephemeris was provided by the Naval Astronautics Group directly and not 
derived from transmitted ephemeris in coded form. Comparisons have been 
made in the position and velocity components in the earth-fixed Cartesian 
system as well as in the spherical system in a procedure similar to that 
described in Section 3.5.1. 

In the case of Data Sets D and S, the precise ephemeris was in the 
NWL 9D system, while the broadcast ephemeris was in the WGS 72 system. 

So a coordinate transformation was performed to bring the precise ephemeris 
into the WGS 72 system before carrying out the comparisons. With the param- 
eters mentioned in Section 3. 3, the equation for the transformation is given by 
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-X- 


■X' 


“AL 

w 

o - 


-x- 

Y 

= 

Y 


-W 

AL 

0 


Y 

-Z- 


-Z- 


. 0 

0 

AL- 


-Z- 


WGS72 9D 9D 

Since the transformation parameters are time independent 



where 

AL = -0.8263 ppm 

W = -0V26 [Vincenty, 1976] 

The broadcast ephemeris was derived from majority-voted ephemeris 
message from data collected at three stations in the case of Data Set D and 
from the coded ephemeris message provided on tape by NAG in the case of 
Data Set S. Where data pertains to an overlap between successive orbit fits, 
the more recent data has been used. The results of the comparisons are shown 
in Tables 3-10 to 3-13. 

Based on the above finding, the following observations can be made about 
the accuracy estimates of the broadcast ephemeris in the WGS 72 system. 

(i) Taking conservative estimates, the positional accuracy of the broad- 
cast ephemeris, taking the precise ephemeris as the standard, can be taken as 
10 m, 10 m, and 17 m in the X,Y,Z directions, respectively. The esti- 
mates for the in-track, out -of plane, and radial directions are 19 m, 15 m, 
and 9 m, respectively. The values obtained in a comparable procedure by 
White [1975], before upgrading of the broadcast ephemeris computational system, 
were 25 m, 15 m, and 10 m, respectively. The smaller values now obtained 
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Table 3-11 


! 

Accuracy Estimates of Broadcast Ephemeris for Velocity Components in Cartesian System 









































Table 3-12 


Accuracy Estimates of Broadcast Ephemeris for Position Components in Polar System 


1 





r iml 1 

mm 

RM S 

SI) 

A * 

RMS 

SD 


RMS 

SD 

■ 

13 


- . 198 

. 625 

. 59:1 

. 04 4 

. 461 

. 450 

sm 


5 . 4 

H 

19 


mm 

. 667 

. 6 12 

. 022 

. 4 54 

. 454 

mm 


6 1 


19 

in 

-.014 

. 14 7 

.14 7 

- . 148 

. 324 

.28 9 

n 

MM 

3. 7 



B9 

. 070 

.20 5 

. 19 3 

. 170 

. 509 

. 482 

mm 

til 

4. 2 


mm 

mi 

. 282 

B9I 

.282 

MM 

. 386 

. 37 1 

>3.3 

5 . 5 

4 . 3 

s 

i 

H 

19 

. 298 

n 

.22 1 

ml 

. 301 

.28 3 

-7. 4 

6 . 3 

3 8 


Table 3-13 

Accuracy Estimates of Broadcast Ephemeris for Velocity Components in Polar System 


n*i& 1 

SaMlIte 

No. of 

1 





filO' 

, 2 

■ (m/s) 

Set 


IMInts 


r- 

RMS 

SD —i 

*tb 

MM S 

SD 

*r 

B*is 

SD 


13 

■ 

496 

.039 

. 258 

.so 

- . 008 

. 028 

.027 

- . 002 

.011 

.01 1 

1. 

19 

506 

- .058 

. 34 9 

.345 

.014 

.03 1 

.028 

. 00 1 

011 

.011 


19 I 

199 

mm 

.13 1 

.13 1 

• 022 

. .380 

MM 

. 00 1 



1) 

1 20 

13 1 

mm 

. 127 ' 

1 .126 

.006 

. 064 

■9 

. 001 




12 

2005 

-03< 

. 280 j 

.278 

| - .016 

. 207 

.207 

- 009 


mmm 

s | 

19 

1879 

.033 

. *i7 4 j 

. 272 

- . 022 

. 286 

. 285 

- . 004 

. 14 8 

. 148 



















































would represent the effect of the improvements incorporated in the ephemerides 
since December, 1975. Previous study by Wells [1974] obtained 26 m, 10 m, 
and 5 m for the above estimates, in a slightly different procedure. 

(ii) The values for the accuracy of the velocity components would be 

0.1 ,0.2, and 0.2 m/sec in the X,Y,Z directions, respectively. The corres- 
ponding values in the in-track, out-of-plane, and radial directions would be 
0. 2 , 0. 1 , and 0. 2 m/sec, respectively. 

(iii) The RMS values for velocity differences are much smaller in die 
case of Data Set L though this refers to a period before incorporation of im- 
provements in the computation of broadcast ephemerides . This is because in 
Data Set L the predicted state vectors of the satellites (both position and 
velocity components) were available directly on cards and referred to a pre- 
coded and preinjection stage, while in Data Sets D and S these have been 
derived from coded message and a fitting process for obtaining time deriva- 
tives of variable elements. So the velocities in Data Set L are free of trun- 
cation errors and errors of fitting. 

(iv) At first sight it may appear irrational that the precision estimates 
are of the same order and, in some cases, even lower than die accuracy esti- 
mates. However, as pointed out in the previous section, the precision estimates 
in reality are the estimates of the prediction errors accumulated over the inter- 
injection period. The overall uncertainty of the broadcast ephemeris as obtained 
by a user will depend on the relative position of the observation time in the 
interinjection period as well as the interinjection period itself, which was found 
to vary from about 9| to 13 1 hours in the data available in this study. 

Thus, if a user happens to track a satellite just before injection, the 
uncertainty of the ephemeris obtainable by him would be a result of the com- 
pounding of the estimates of precision and accuracy given above. In this case, 
the positional uncertainty of the satellite is likely to be 26 m in-track, 20 m 
cross-track and 10 m in radial directions. However, for a user tracking a 
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satellite during or soon after injection, the uncertainties will be those given by 
the accuracy estimates above (i .e. , about 19m in-track, 15 m cross-track, and 
9 m in radial directions). If prediction errors are not taken into account, the 
improvements in the Transit system in December, 1975, have brought the broad- 
cast ephemeris closer to the precise ephemeris. 

(v) Another observation which can be made is about the prominent nega- 
tive value of Ar (between three to seven meters) in all data sets. 

The possibility of a radial bias is evident and has been pursued in 
Chapters 4 and 5. 

(vi) The existence of cross-track and in-track biases cannot be complete- 
ly overruled. As regards cross-track bias. Data Set D indicated a bias of about 
-0701 for satellite 19 and 0'.'07 for satellite 20. In Data Set S both satellites 12 
and 19 indicated a bias of about 0728 and 0730, respective^ . While in Data Set 

L the bias indicated is of the order -0720 and -0728. A similar situation is seen 
in the case of in-track bias . This led to the conclusion that in the data sets used 
in this study, a consistent evidence of the existence of in -track and cross-track 
bias is not available. These biases have, therefore, not been pursued further 
in this study, but a need does exist for identifying these biases with more data 
sets in the future. 

This conclusion was also supported by the results of tests of hypothesis. 
Only in the case of radial bias the tests indicated that the hypothesis that the 
expectation of the bias is zero could be rejected at a =0.05, in the case of all 
data sets. The biases investigated here are in the nature of constant systematic 
effects that can be associated with the broadcast ephemeris irrespective of the 
satellite or its pass tracked. Other biases in individual satellite passes, if 
varying in magnitude and sign between passes, may average out to an insig- 
nificant value over a large number of satellite passes . These latter types of 
pass biases, referred to in the literature [Kouba and Wells, 1976], are viewed 
as signals which can be separated in an adaptive filtering procedure discussed 
in Chapier 5. 
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Based on the above, and taking into consideration the 2m uncertainty 
in each component of position expected of the precise ephemeris, it can be 
stated that the positional uncertainty of ephemeris as obtainable by a user just 
before injection is about 35 m. The positional uncertainty obtainable immed- 
iately after injection is about 25 m. The figures given in the literature places 
the single pass error budget for a surveyor at 12-28 m for the broadcast 
ephemeris in WGS 72. 
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4 . STATISTICAL ANALYSIS O F EPHEMERIS DATA 


4.1 Introductory Remarks 

An immediate application for the accuracy estimates of the broadcast 
ephemeris, obtained in the study described in Chapter 3, lies in utilizing this 
information in the form of a diagonal matrix to assign appropriate weights to 
the orbit parameters in an adjustment procedure in which the main parameters 
of interest would still be the station positions. 

In this method it is assumed that the uncertainty in the orbit is the result 
of a system "noise" arising from a random variable whose outcomes are inde- 
pendent and identically distributed with a zero expectation and finite variance. 
This is similar to the assumption invariably made about the observational 
"noise" or random observational "errors" in any least squares adjustment pro- 
cedure, irrespective of whether a system "noise" is taken into consideration or 
not. In effect, the procedure amounts to adding the best known orbit parameter 
values as additional observations in the adjustment. 

Though the mathematical treatment for system "noise" and observational 
"noise" is similar in the above situation, the distinction between the two is clear. 
A system "noise" arises from an inadequacy in the mathematical model, while an 
observational "noise" arises from the measurement process. Further in this 
study, system "noise" arises from two models, the model which describes the 
dynamics of satellite motion and the model which relates the observable quantity 
with the parameters. The first model is the concern of this chapter. 

In many situations, such as the satellite orbit in this study, the succes- 
sive outcomes of system noise are correlated, though their expectation may 
still be zero. In such situations estimates of the parameters of interest may 
improve if the system "noise" is modeled by a stochastic process in the form 
of a "signal" to be estimated along with ihe other parameters, instead of being 

45 



compensated for by weighting parameters . 

The situation is well described by considering a physical system whose 
dynamic behavior can be modeled by the first-order linear differential equations 

$ = A(t)5f + G(t)W(t) fort ^tb where X = X(t) (4.1) 


Here X is an n vector called the state of the system. We say that X is 
a state vector if X (tj) can be determined unambiguously from a knowledge of 

X (to), ti s t 0 , and A(t), G(t), and W(t) for tb 5 t s ti. X contains all the 
parameters of interest required to describe the system. W is an n x 1 

"disturbance” vector, a term which implies a "noise" if the outcomes are 
uncorrelated and a "sigr/il" if correlated, and t denotes time. A (t) and G(t) 
are n x n system matrices assumed to be continuous in time t. It is assumed 
that the initial time t 0 is fixed and initial state X (t Q ) is known. The forms of 
the matrices A(t) and G(t) are also given. 

The observation process is modeled as 

Z(t) = H (t) X (t) + V(t) (4.2) 


where Z is a p vector of measurements , V is a P vector of measurement 
errors, and H(t) is a continuous p x n design matrix or measurement matrix 
of known form. A block diagram for the process implied by equations (4.1) 
and (4.2) is shown in Fig. 4.1 as adapted from [Meditch, 1969]. 



Fig. 4.1 Block diagram for continuous linear system description. 


46 


ORIGINAL PAGE'S 
OF POOR QUAUm 




Estimation procedures in the above system require some knowledge about 
the statistics of the disturbance as we 11 as the observational noise. In most cases 
where the former is not directly observable, observational residuals are the 
only source of information for determining the statistics of the disturbance as 
well as of the observational errors, and it is difficult to separate their effects. 
However, there are situations where disturbance is accessible to form some esti- 
mate about its statistics. Availability of.precise and broadcast ephemerides of 
the same satellite gives rise to one such situation. The outcomes of this distur- 
bance process are generally correlated, and it can more appropriately be term- 
ed as a signal. 

Precise and broadcast ephemerides provide us information about the state 
of the same satellites at two-minute intervals for observation spans over a time 
period. Precise ephemeris is known to be much more accurate than broadcast 
ephemeris. So treating precise ephemeris as the standard, t»u small differences 
obtainable by comparing the two can be treated as the outcomes of the system 
signal inherent in broadcast ephemeris and their statistics estimated by the nor- 
mal sampling procedures. Immediately, several possibilities arise as these 
comparisons can be carried out in different ways. The specific procedure for 
comparison depends on how the state disturbance is proposed to be viewed. For 
example, the state disturbance could be viewed as predominantly a consequence 
of uncompensated acceleration, in which case a stochastic process would be re- 
quired to be added to the relevant first-order differential equation for accelera- 
tion in the dynamic model. Outcomes of this process would be computed by 
comparing the accelerations given by precise and broadcast ephemerides of the 
same satellite, at the same epochs, and in the same coordinate system for 
obtaining their statistics. Further, the uncompensated acceleration could be 
considered either in the Cartesian system or a spherical system. 

On the other hand, the state disturbance could be viewed as predominant- 
ly a consequence of a positional or a velocity signal associated with broadcast 
ephemeris, which implies that the uncompensated accelerations are negligible, 
but due to some reason, there is a systematic deviation in position or velocity 
which can be associated with the broadcast ephemeris as a stochastic process. 
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In this investigation four different types of signals have been studied 
by computing their statistics, and one type of signal has actually been used to 
study its effect in local orbit refinement. 

The theoretical stochastic process of the signal is denoted as £(t), while 
the practical outcomes of this process are denoted as s(t) where t denotes time. 
The autocovariance of £(t) is denoted as R(u), while the sample autocovariance 
obtained from data is denoted as C(u), where the lag parameter u represents 
the separation in t. The subscripts clarify the specific signal intended. The 
signals and their statistics, as obtained in this study, are described in this 
chapter, and their applications are discussed in Chapter 5. 

4.2 Computation of the Signal Outcomes 

The ephemerides provide position and velocity components in the earth- 
fixed systems. As indicated in the previous chapter, these components are 
more correlated than the in -track, cross-track, and radial compone nts 
which are obtainable by transforming the state vectors into a spherical system. 
As explained in Section 4.1, keeping precise exhemeris as the standard, reali- 
zations of the following four signals were computed for Data Sets D and S which 
are in the current WGS 72 system: 

(i) Acceleration signal £^(t), £y (t), £^(t) in the Cartesian system 

(ii) Acceleration signal £^(t), £^(t), £j:(t) in the polar system 

(iii) Velocity signal £^(t), £ r - (*) in the polar system 

(iv) Position signal 4^(t), ^r ^ ^ the polar system 

All signal outcomes were computed by first transforming the precise 
ephemeris state vectors to the WGS 72 system of the broadcast ephemeris. 
However, it has been found both analytically and numerically that the transfor- 
mation has a negligible effect on the acceleration signal. This is shown in 
Appendix A. 

The actual mechanics for the realizations of the above signals will 
now be described, followed by the procedure used to obtain the statistics. 
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4. 2. 1 Acceleration Signal 


If Xp(ti) and X p (ti +2 ) are the X components of a satellite velocity in 
m/s according to the precise ephemeris after transformation to the WGS 
72 system, and X„(ti) andX n (ti +2 ) are their corresponding values according 
to the broadcast ephemeris at even-minute epochs ti and ti+a, then 


Xp(t, +1 ) 


and 


X # (ti + i) 


Xp(t, +2 ) - X p (ti) 
120 


X a (t 1+2 ) - X a (t t ) 
120 


in m/s a 


in m/s 3 


(4.3) 


(4.4) 


The realization s- (ti +1 ) of the acceleration signal 4 - (ti+i) will be given by 
X X 

s~ (ti+i) = X p (tj+i) - X a (ti+i) (4.5) 

A 

Similar computations yield the realizations s^. (t) and s- (t) of the signals 
4y(t) and £z<t). 

To obtain the outcomes of the signal 4 x (t), 4^(t)* 4f(t)» a similar 
procedure was used by first transforming the precise and broadcast ephemeris 
state vectors from the Cartesian system to the spherical system as explained 
in Section 3. 5. 1. 


4.2.2 Velocity Signal 

Similarly, if A p (ti) and X a (ti) are the components of a satellite velocity in 
"/s according to the precise and broadcast ephemeris, respectively, at the 
same even-minute epoch (ti), the outcome (ti) of the velocity signal 4 X (ti) 
will be given by 

s x (t») = Xp(ti) - X a (ti) (4.6) 

Similar computations yielded 

s^(tt), s^(ti) the realizations of the velocity signals 
4f.(tt) in the spherical system. 

49 „ i 

ORIGINAL PAGE 

OF POOR QUAtfl 


4. 2. 3 Position Signals 

Following the same procedure, the realizations s^(ti) of the signal 
were obtained from the relation 

s x (ti) = X p (t t ) - X a (ti) 

where X p (ti) and X n (ti) are the X position components of the satellite at oven- 
minute epoch ti in the precise and broadcast ephemeris, respectively. Similar- 
ly, the realizations s<p(tj), s r (ti) of the signals £<p (ti) and £ r (ti) were obtained. 

4.3 Computation of the Signal Statistics 

Considering each of the signals £(t) above as a stochastic process, the 
computations of Section 4.2 yield the discrete observations s(ti) or outcomes of 
the continuous time series at equal intervals of two minutes over the available 
data spans. 

For computing further statistics, the following assumptions were made: 

(i) The process £(t) has reached a steady state in the sense that the 
statistical properties of the series are independent of absolute time. This 
implies that the probability density function of each signal is independent of 
time, has a constant mean fiand constant variance a 3 . That is, it forms a 
stationary time series. 

(ii) The process exhibits the property of ergodicity, which enables the 
computation of a time average over a record to represent an ensemble average. 

With these assumptions the outcomes s(t) can be used to form the sample 
autocovariance function, the sample cross -covariance function, and the sample 
autocorrelation functions. 

For example, considering £ " (t), its autocovariance function is 

A 

defined as 

Ii££(ti, fc) E[($a<ti) - pgft)) (ta) - M*‘(ta) ) ] (4.6) 
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where E is the expectation operator, and (t) = E (t)] . Because of the 
stationarity assumption, R^(ti, tg) immediately becomes a function of the 
time lag u = tg - ti only, and 

( u) = E - M*) (€*(t+u) - fi*)] (4.7) 

Sample autocovariance C~ ~(u), its estimate from data is computed from 
the outcomes s»(tj) using the relation 

X 


N-U 

C XX^ = " s x> f®X ^ tl+u ^ ' 


where 


** “ n5 


i represents the epoch of each data point taken at two- minute 


intervals in an ordered set 
N is the total number of data points 


[Jenkins and Watts, 1968] 


The sample autocorrelation function r^(u) is given by 


.. c xx w 

r XX U “ Cjtx(o) 


Similarly, if we consider two signals £ - (t) and £ - (t), the cross -covariance 

X Y 

function for lag u becomes 

R*Y < u > = E (t) ‘ M X> ^ y (t+ U) “ M Y } 1 

for which the sample cross -covariance function is 

N-U 

C XY (U) = nE (S X (tl) ' J x )(S Y (tU)i) ~*Y> (4 ‘ 9) 

i= l 

where 


ijE S Y (, ’> 
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For equations (4.8) and (4.9), an alternate formulation is also available 
in the literature [Anderson, 1971] according to which the divisor (N-u) is used in- 
stead of N for computing the estimates. Divisor N has been used in this investi- 
gation as estimators with divisor N usually have smaller mean square error and 
are positive definite [Jenkins and Watts, 1968 p. 184]. For lag o, the sample 
cross-correlation matrix for signals 4 " (t), 4y (t), 4” (*) is formed by unities 
along the diagonal and factors like 



(o) 



Cyy (o) ) 


as off-diagonal terms. 


Similar statistics for all signals have been obtained and are shown in various 
figures as follows: 


Fig. 

4.2 

c xx <u> 

for Data Sets 

D and S 

Fig. 

4.3 

r xx (u) 

for Data Sets 

D and S 

Fig. 

4.4 

C YY (U > 

for Data Sets 

D and S 

Fig. 

4.5 

r~ - (u) 
YY ' ' 

for Data Sets 

D and S 

Fig. 

4.6 

c zz <u) 

for Data Sets 

D and S 

Fig. 

4.7 

r zz (u) 

for Data Sets 

D and S 

Fig. 

4.8 

C XA (U > 

for Data Sets 

D and S 

Fig. 

4.9 

’U (u > 

for Data Sets 

D and S 

Fig. 

4.10 

c „<p <“> 

for Data Sets 

D and S 

Fig. 

4.11 

r . (u) 
< p<p ' ' 

for Data Sets 

D and S 

Fig. 

4.12 

C (u) 

for Data Sets 

D and S 



rr 



Fig. 

4.13 

r rr< U > 

for Data Sets 

D and S 

Fig. 

4.14 

C tt <“> 

for Data Sets 

D and S 

Fig. 

4.15 

r U < u) 

for Data Sets 

D and S 

Fig. 

1.16 

C> • (u) 
<P<P 

for Data Sets 

D and S 

Fig. 

1. 17 

r. . (u) 
<p<p 

for Data Sets 

1) and S 

Fig. 

1. 18 

C..(u) 

rr 

for Data Sets 

D and S 

Fig. 

4.19 

r. . (u) 
rr 

for Data Sets 

Dand S 
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U (minutes) 


Figure 4.2 Sample Autocovariance C^^(u) 
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Figure 4.3 Sample Autocorrelation r^^(u) 
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Figure 4. 4 Sample Autocovariance C^(u) Figure 4.5 Sample Autocorrelation r^(u) 
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Figure 4 . 10 Sample Autocovari ance C (u) 
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Figure 4. 11 Sample Autocorrelation (u) 
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. 14 Sample Autocovariance Ci i (u) Figure 4. 13 Sample Autocorrelation r;j ^ u ) 


















Fig. 4.20 

°u (a) 

for Data Sets D and S 

Fig. 4.21 

r u<"> 

for Data Sets D and S 

Fig. 4.22 

c w (u) 

for Data Sets D and S 

Fig. 4.23 

r~ - (u) 
<P<P 

for Data Sets D and S 

Fig. 4.24 

C--(u) 

for Data Sets D and S 


rr 


Fig. 4.25 

r rr^ 

for Data Sets D and S 


The sample cross-correlation matrix for the signals has been worked out for 
a typical case of satellite 19 for the lag u = 0 for Data Set S. The matrices 
obtained are as follows: 

Sample Cross -Correlation Matrix for £-(0), £-(0), £-(0) 

Satellite 19 X Y Z 

V X S Y S Z 

'1 0.589 -0.073' 

1 -0.234 

1 


Sample Cross -Correlation Matrix for £. (0), £ (0), £ (0) 
Satellite 19 r 


n 


^0 

-0.018 


^ r 

-0.1681 


1 


-0.135 

1 


Sample Cross -Correlation Matrix for £:(0), £ • (0), £j.(0) 
Satellite 19 ^ 


ri 


- 0.011 


’ r 

0.029 


1 


-0.003 

1 
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Figure 4.20 Sample Autocovariance c a <u » 


Figure 4.21 Sample Autocorrelation r££(u) 






























Sample Cross -Correlation Matrix for £c-(0), £ -(0), £ --(0) 
Satellite 19 A <p 



0.041 

1 


*r 

0.005- 

-0.042 

1 


Table 4-1 gives a summary of the data analyzed in obtaining the 
statistics . It may be noticed that satellite 19 is common to Data Sets D and 
S. The number of data points shown refer to the number of realizations of 
the acceleration signal. The realizations for other signals are larger than 
this number. 


Table 4-1 

Summary of Data Analyzed 
for Signal Statistics 


Data Set 

Satellite 

Minimum No. of 
Data Points 

D 

19 

182 


20 

118 

s 

12 

1674 


19 

1502 


4.4 Observations 

The purpose of obtaining sample autocovariance functions and related 
quantities is to determine a model for the underlying stochastic process, 
in this case, the various signals. The following observations can be made 
based on the values obtained taking the various types of signals in tui-n. 
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4,4.1 Acceleration Signals (u), £^.(u), ££( u ) 

(i) The signal with the smallest autocovariance has a magnitude of over 
2.0 x 10” 8 m 3 /s 4 . Acceleration of 1.0 x l(f 4 m/s 3 can cause a positional devia- 
tion of 1.8 m after ten minutes which can be taken as the observation period of 

a pass. So the signal is not insignificant where submeter accuracies are sought. 
The largest covariance is over 1.0 x 10”* m 3 /s 4 . 

(ii) Except for satellite 12 of Data Set S, the autocovariances are posi- 
tive for a lag up to four minutes. For £ - (u) and £ ~ (u) of satellite 12, the auto- 

Y » 

covariance is positive up to a lag of two minutes . 

When dealing with Doppler observations with a JMR-1 receiver as in 
Data Set D, the Doppler counts are available at every 4. 6 seconds approxi- 
mately, and if three stations coobserve a pass, data is acquired at a rapid 
pace of about 78 observations for every two minutes observed in common. So 
this signal is worth attempting to model in an effort to refine the orbit. 

It will be seen in Chapter 5 that a simple and mathematically tractable 
model is sought for the stochastic process. A simple and convenient model for 
the autocovariance is an exponential form without any periodic term whieh gives 
rise only to nonnegative values. The lag up to which the sample autocovariance 
is positive is therefore very relevant. 

(iii) The sample cross -correlation matrix of component signals indicates 
correlation of about -0. f,3 between £ r, and £ s and a correlation of about 0.59 

1 /j 

between £ - and £ -. So this is not the ideal set of signals to attempt to model 

A Y 

in any procedure where an assumption of noncorrelation of the component 
signals simplifies covariance propagation. However, it is not too high either 
to leave this procedure completely out of consideration. 


67 




4.4.2 Acceleration Signals £ ^(u), 4 -( u ). 4f( u ) 

(i) The signal with the smallest autocovarianee has a magnitude of 
over 0. 2 x 10“ 10 ("/s") which can cause a positional deviation of one meter 
after ten minutes. The largest autocovariance has a magnitude of over 1. 2 
x 10 -6 m/s 4 . 

(ii) Regarding the lags up to which the autocovariance is positive, the 
same remarks as for £ -;(u), £ -(u), £ a(u) above apply in the case of £ - (u) and 

A Y I* CP 

£ - (u). For the signal 4c( u )t the positive autocovariance persists even after 
F A 

eight minutes. 

(iii) The sample cross-correlation matrix shows the largest correlation 
as low as 0. 01 making this set of signals a good candidate to model. 

4.4.3 Velocity Signals £ ^ (u), 4^ (u). 4f( u ) 

(i) The signal with the smallest autocovariance is over 0. 4 x 10 _B ("/s f. 
A velocity of 0. 1 x l(f 3 "/s implies a positional deviation of over 2 m after 

ten minutes. The maximum autocovariance is in the order of 14. 0 x 10” 6 ("/s) 2 . 

(ii) Regarding the duration of positive covariance, the same remarks as 

above apply for signals £ • (u) and £ . (u). In the case of the signal £ : (u), the 

0 r A 

positive correlation continues to six minutes. 

(iii) The largest correlation is only 0.03 as seen from the sample cross- 
correlation matrix making this another good candidate for an attempt to model 
the signal. 

4.4.4 Position Signals £ . (u), £ (u), $ (u) 

A <P 1 

(i) In this case the smallest autocovariance is over 0.02 (") 2 . O'.'l 
implies a deviation of 3. 6 m at the altitude of the satellite. The largest auto- 
covariance is over 0. 2 ("f. 
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(ii) The duration of positive correlation is much larger than in the cases 
of all previous signals. It persists at least up to eight minuter in case £ r (u) 
for satellite 19 in Data Set S and extends over 18 minutes in other cases, making 
this a good candidate for an effort to model. 

(iii) The maximum correlation indicated in the cross-correlation matrix 
is 0.17 which is fairly low making this signal worthwhile being modeled. 

(iv) Another interesting observation in the case of these auto covariances 
is the value of the mean of the signal realizations obtained in forming sample 
autocovariances C rr (u). The mean in all cases was a negative value varying in 
magnitude from -3.3 to -7. 8 m. The effect of treating this as a bias in a point 
positioning experiment has been described in Chapter 5. 

4.5 Application in Present Study 

The computations of the realizations of the signals and their sample 
statistics have been described above. Some general observations about the 
statistics obtained have also been made. 

It can be noticed that satellite 19 is common to Data Sets D and S, and 
Data Set S is over five times larger than Data Set D. So satellite 19 was a good 
candidate for experimentation, giving a greater weight to the statistics obtained 
from Data Set S . 

In each of the figures 1 .2 to 4.25 there are two curves witli respect to 
satellite 19, one pertaining to the period January, 197C>, and the other to the 
period October, 1976. It is comforting to see that these two curves are close 
to each other in most figures indicating a nearly stabilized situation. 

It will be seen in Chapter 5 that the model chosen for the stochastic 
process of the signal £(t) is of the form 

i(t) -0 4(t) + o0 w(t) 
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where 


w(t) is the Gaussian noise with the properties Efw(t)] -- 0, 

E[w(t) w(s)] - 5(t - s) 
a is the variance parameter 

6 (t- s) is the Dirac delta fune' ion 
j8 is the time correlation coefficient 

This choice arises from the fact that if it is further assumed that £ (o) ~ 

N(o, a 2 j3/2), independent of lw(t)} in the above model, then 

R^(u) = a 2 (0/2) e“^ (U) [Jazwinsky, 1970] 

which is the analog of Uie first-order autoregressive process or first-order 
Gauss-Markov process. In this situation, the curves obtained in this study 
can directly give the initial estimates of the parameters of the autocovariance 
function in procedures to be described in Chapter 5. 

It is evident that with the material available many different experiments 
and approaches are possible, depending on which signal is proposed to be 
chosen to represent the state disturbance. In this study, however, the con- 
straints on time available to the author and the form of software available have 
made only some types of experiments possible. These have been discussed in 
the next chapter. 
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5. EXPERIMENTS FOR IMPROVED POSITIONING 


5. 1 Introductory Remark s 

With the statistical information obtained in Chapters 3 and 4, a couple 
of experiments were undertaken in which this information was utilized. These 
experiments and their results have been described in this chapter. 

At this stage it is appropriate to indicate the overall philosophy adopted 
in designing these experiments, for the ultimate aim of gaining improved 
station positions with broadcast ephemeris and Doppler data. 

The Doppler observational data available and used pertains to Data Set 
D in which three stations have coobserved several passes . As discussed in 
Chapters 2 and 4, it is possible to carry out a solution in a short arc mode 
and obtain values for both the orbit unknowns and the station unknowns. But 
as indicated in Chapter 2, such solutions may have a problem which is pre- 
cisely stated by Mueller [1976]: 

The results are being scrutinized by theoreticians who regard the 
results as "meaningless" in view of the fact that the dynamic solu- 
tions are rank deficient and as such the problem (of simultaneously 
determining geocentric station and satellite parameters) is theo- 
retically unsolvable, e.g. , the system of reference defined by such 
solutions would depend entirely on the a priority selected values of 
parameters (e.g., station coordinates). 

In this study, the knowledge of the satellite positions is available in the 
form of broadcast ephemeris. The broadcast ephemeris gives a set of values 
in a consistent coordinate system with uncertainties which have been assessed 
in the study and described in Chapter 3. 

Doppler observational data from three coobserving stations (Data Set D) 
has been used with broadcast ephemeris to obtain the best values for the 
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station coordinates which eould bo held fixed in a subsequent filtering program 
for orbit refinement. Under this situation, and keeping in view the rank 
deficiency problem described earlier, a decision was to be taken about the 
specific constraints to be used in this study in the solution for station positions. 

As described by Pope [1971], the rank deficiency of the normal matrix 
can be obviated by use of either (a) weighted constraints or (b) absolute con- 
straints on parameters. Method (a) becomes (b) as weights on the parameters 
are increased to infinity. 

As regards the number of constraints, they can be either the minimal 
set of constraints, or they can be more than minimal. Among minimal con- 
straints, there is the attractive possibility of using an inner constraint solu- 
tion |Blaha, 1971] which leads to the "free network adjustment" without 
explicitly forming a pseudo- inverse matrix. 

Although inner constraints yield an optimum coordinate system with a 
minimum trace for the variance-covariance matrix of the parameter esti- 
mates, the solution and the coordinate system itself is optimized with respect 
to the initial values of the parameters used in the adjustment. While this may 
be advantageous in a larger network adjustment where the coordinate system 
also needs to be optimized, in a local adjustment like the one used in this 
study where a coordinate system in unambiguously accessible (viz. , the \V(’>S 
72 system through the broadease ephemeris), it is best to adopt it and to 
ensure that the values obtained in the adjustment are in the chosen coordinate 
system. 

Since the ground station positions are the main unknowns, this could be 
achieved bv absolutely constraining the satellite orbits by introducing very 
large weights. Thus it was decided to hold the satellite orbit fixed in the solu- 
tion for station |K>sitions with SAHA. 

Adopting and maintaining the WHS 72 coordinate system as described 
above, the problem is tackled in the following three stages: 


(i) Determination of station positions with observational data, holding the 
orbits fixed, 

(ii) Improving the orbits keeping the station positions unchanged, 

(iii) Redetermination of station positions holding the improved orbits fixed. 

Within the time available, only two experiments could be completed, 
and only these two experiments have been described. In the first experiment 
the radial bias noticed vide Section 3.5 has been taken intoaccount, and the effect 
of removing this suspected bias in the broadcast ephemeris on station position 
recovery has been determined. This experiment was aimed only to study the 
radial bias identified in Chapter 3. Similar studies must be carried out with 
respect to other biases, but for reasons explained in Section 3.2, these have 
not been pursued in this investigation. As a follow up of this experiment, the 
reason for the radial bias has also been investigated. 

In the second experiment, the first steps have been taken to judge the 
feasibility of utilizing the statistical information obtained in Chapter 4 in an 
adaptive filtering procedure for orbit improvement. In this experiment, the 
station positions have been held fixed at values obtained in stage (i). Statistical 
information obtained earlier was utilized, but simulated range rate observations 
available from existing software were utilized to judge the performance of the 
filtering procedure. 

Simulated range rate observations have been used in this study only as 
a first step to judge the feasibility of using the adaptive filtering approach in 
practice, so that the available software could be used with minimum modifica- 
tions. In position determination with Doppler data, range differences would 
have been the appropriate observations. 

Range difference observations have a more direct influence than range 
rates on the position of a satellite. So in using simulated range rates to judge 
the feasibility of the filtering pi'ocedure, an observational mode less optimum 
than that obtained in reality is being adopted. The next steps in this study 
would have been to first try the procedure with simulated range differences 
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and then to use real data . It is best to undertake these steps after first 
deciding the optimum model for the state disturbance. For example, the 
formulations for viewing the state disturbance as an acceleration signal in the 
Cartesian coordinate system would differ considerably from the formulation 
for state disturbance viewed as a position signal in the spherical coordinate 
system. Unfortunately, in the time frame available to the author, it was not 
passible to complete these further steps. 

The presence of state disturbance and its role in the evolution of the 
shite characterizes the main distinguishing feature of the filtering procedure 
in this study. But the method of estimation from observations is similar to 
the least squares techniques. Thus it was felt appropriate to review the 
mathematical formulations leading to the conventional sequential filter 
algorithm ( also called the first order filtering technique! as a back- 
ground to the adaptive filtering procedure in the second experiment. 

These are included in Appendix B. 

f> . 2 Kxpcriment for Radial Bias 

5 .2.1 Detenu ination of Stnt ion Posit ions 

Station positions for three stations in the Morula area were obtained with 
the help of the Short Are Geodetic Program (SAGA) (Brown :uid Trotter, H><*9) 
available at The Ohio State University. The formulation used in the program is 
documented both in the above- referred publication (Kumar, 197(1] and in (Brown, 
197<>] . 

With the version of SAGA available at OSU, the following main steps 
were required to obtain a solution: 

(i) modification of the available subroutine for Geoceiver raw data to make it 
compatible with JMR raw data, 

(ii) majority vote of the ephemeris message to obtain state vectors of 
satellites at two-minute intervals, 

(iii) computation of mid arc state vectors for each pass for input to SAGA 
with the program SAMVAP. 
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The meteorological data acquired at the stations formed part of the input. 

The solution was obtained from 12 passes holding the orbits fixed and 
allowing the station positions to remain free in the adjustment. 

5.2.2 Determination of Radial Bias 

As mentioned in Section 4.5, a negative radial bias was noticed while 
computing the sample autocovariance C rr ( u ) for both Data Sets S and D. To 
assess a more reliable value of the suspected bias, Data Sets S and D for the 
four satellites were combined with the following results: 

degrees of freedom 4313 
bias = A r = -5.3 m 

Otr - 4.5 m 

A t-test is not strictly valid because of the fact that adjacent outcomes 
of A r are correlated. However, in view of the large degree of freedom avail- 
able, a t-test w.'S carried out to test the hypothesis Ho: M = 0 against Hi: fi 4 0 
at or = 0.05. 

-5 257 

t (computed) = 4 . 50 7//4Sl4 ' - 76 ' 61 

t (tabular) = 1.960 

Therefore the hypothesis that fi = 0 is rejected. 

5. 2. 3 Correction of the State Vectors for Bias 

The mid arc state vectors used in Section 5.2. 1 were corrected for a 
radial bias of -5.3 m. As shown in Section 4.2. 3, the radial signal s r - 

r p - r n = Ar. The bias is the mean value of s r . The corrected radial distance 

r c = r n + Ar. 

For each pass used in the solution in 5.2. 1, the mid arc state vector 
was first transformed from the Cartesian components (X,Y,Z) to the polar 
components (X, O, r). The radial component was corrected to r c . The corrected 
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components (X , <p , r c ) were transformed back to the Cartesian components (X c , 
Y c , Z c ), and the solution of Section 5.2.1 was repeated. 

5. 2. 1 {^determination of Station Positions After Removal of Bias 

Table 5-1 gives the results of station position determination obtained 
before and after removal of radial bias. Table 5-2 shows the weight coeffi- 
cient matrix which remained unchanged to the number of digits shown. 

5.2.5 O bs ervati ons 

(i) The a posteriori variance of unit weight was reduced by a very small 
amount (.0003) after removal of bias. The weight coefficient matrix is unchanged 
to four decimal places. The station coordinates changed by very small amounts 
(< 0.1 m). 

(ii) The bias was absorbed by small changes in the values of statical coor- 
dinates and other pass parameters such as frequency bias, frequency drift, and 
refraction scaling factor. 

(iii) One external check was available. The distance between stations 1 
and 3, according to terrestrial survey, was 29 360.880 m, as given in (Brown, 
1976]. In the above determinations, the distance was 29 361.003 m before cor- 
rection for bias and 29 360.928 m after removal of bias. 

The large standard deviation of the chord recovered (8.7 m) pre- 
cludes a more pos itive statement . However, the removal of bias has made 
the solution for chord length closer to the terrestrial value. 

5.2.6 Explanation for Radial Bus 

As a follow up of the above experiment, an effort was made to understand 
the reason for the existence of the radial bias by carrying out the following 
adjustment. 

It was assumed that the transformation parameters given by Andcrle 
[1976] between the NVVL 9D of precise ephemeris and the WGS 72 of broadcast 
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Table 5-1 


Results of Station Position Determination 

State Vectors As Broadcast : 

Apriori Standard Deviation ofObservation = err =20 cm 

Degrees of Freedom = 1311 
Final Station Coordinates 


Station 

No. 

X 

(m) 

Y 

(m) 

Z 

(m) 

1 

9 20 731. 202 

-5578835 . 827 

2941252. 526 

3 

892362. 893 

-5579450.679 

2948797.691 

4 

8 85 656. 500 

-5573826. 292 

2961347.463 


& 0 = 1. 07703 


i State Vectors As Corrected for Bias : 

! . 

{ A priori Standard Deviation of Obsei'vations = ar = 20 cm 


Degrees of Freedom = 1311 
Final Station Coordinates 


Station 

No. 

X 

(m) 

Y 

(m) 

Z 

(m) 

1 

920731. Ill 

-5578835.925 

2941252. 572 

3 

892362. 866 

-5579450 . 681 

2948797.692 

4 

885656. 456 

-5573826.278 

2961347.440 


= 1.0767 1 
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Table 5-2 


Weight Coefficient Matrix for Station Coordinates 
( m 2 ) 


Xz 

r.* i *7o+o? 

Yz 

Zz 

x 3 

y 3 

z 3 

x 4 

G.*704n*C? 

C. 76740*02 

— • - 

• - - 

- 

- ... 

. . 

-o.i* n^+c? 

-0.2055040? 

C. 136 1 0402 





3.1 4J40-G6 

0.561 ? D-07 

-0. 16690-05 

C.?4«7D*02 

.... . 

■ — 

— — 

0.1 1073-07 

0.24000-06 

0.51930-06 

-0.95700401 

0.21520402 



-0.47813-07 

0.11600-06 

0.23930— 06 

b.VsoooVo?' 

-0.1 66 704*02 

“ 0.2309040? 

- 

0.23070— 06 

-0. 22670-07 

-0.1 76 30-07 

0.1508r«-06 

-0. 95740-07 

-0.12230-06 

0.21 160*402 

-0.68660-07 

0.2 62 30 - 06 

O. 3^7720- 06 

-0.7849^-07’ 

0.60200-06 

0.3 10? 0**06 

^o75 9ndV6F 


Y 4 Z 4 


6^19*30*02 


0.1Ci5CD-D<> -0.1 2 521-C7 0.353*0-06 -0.33270-07 O.162*0-0ft 0. ^3090-67 0.T 4 82D*O1 -0.10540*02 0^11710*02 


ephemeris take care of the rotations completely but not scale and origin shift. 
With this assumption, a four-parameter transformation was carried out between 
the precise ephemeris, transformed to WGS 72 as per parameters given by 
Anderle [1976] and broadcast ephemeris. Points were selected at intervals of 
over 72 hours to break possible correlation existing between adjacent data 
points. From Data Sets D and S, 26 points were available for satellites 12, 

19 and 20, thus giving 78 observations in an observation equation model. 

Adjustment was carried out with the following mathematical model in 
which AX, AY, AZ and Al were considered as the translation and scale 
parameters 


'X' 


'AX' 


"i + Al 

0 

o - 


-x- 

Y 


Ay 

+ 

0 

1 + AL 

0 


Y 

-Z. 


-AZ- 


- 0 

0 

l + Al. 


-Z_ 


N P 


where 


(5.1) 


N denotes the NAG broadcast ephemeris 

P denotes the precise ephemeris transformed to WGS 72 as per given 
transformation parameters [Anderle, 1976] 


The results obtained were: 


The correlation matrix was ; 


A 

A X 

a 

AY 

A 

AZ 

A 

AL 


0.913 



Ax 

Ay 

AZ 

Al 

-2. 3 

Q > 

r» 

1.4 m 

1 

0. 252 

-0.160 

0. 337 


* 



1 

-0.355 

0. 748 

9.8 

a ; v - 

2.0 m 






AY 




1 

-0. 175 

-2. 0 

a 

q i _ 

2.6 m 






A Z 


_ 



1 

1.43 

> 

fT * - 

A L 

0. 34 ppm 
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A t-test was carried out attt = 0.05 to test if the parameters could be 
considered as 0. The test failed with respect to A Y and AL. So there are 
strong indications that there is an origin shift of about 9.8 m in the Y-axis 
direction and a scale correction of about 1.4 ppm between the two coordinate 
systems considered above. These two parameters can be explained as follows: 

If the origin of the precise ephemeris coordinate system is assumed 
to coincide with the geocenter, the origin of the broadcast ephemeris system 
gets an offset of about 9.8 m in the direction of 90° W longitude, which is 
understandable as all the four stations generating the broadcast ephemeris 
are located in continental United States . 

As regards the scale parameter, it is conjectured that this scale factor 
consists of two components. One component of about 0.9 ppm is due to the scale 
difference always noticed between the Doppler system and terrestrial surveys 
for some unknown reasons. As seen in Chapter 2, the scale in a Doppler system 
is derived from X , the wave length of transmission. The second part of about 
0.5 ppm is the inherent scale difference between the precise and broadcast 
ephemerides . A similar scale factor of about 0. 4 ppm had been detected by 
White [1975] when he carried out the transformations in his study. During that 
time, the coordinate system of both the ephemerides was assumed to be the 
same (NWL 9D). 

5.3 Experiment in Adaptive Filtering 

The aim of this experiment was to assess the feasibility of utilizing the 
statistics collected in Chapter 4 in an adaptive filtering procedure described 
in [Myers, 1973] for orbit improvement. The mathematical formulation of 
adaptive filtering and the description of the software which motivated the present 
experiment will first be given before describing the experiment and the results. 
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5.3.1 Adaptive Filtering 

In the classical sequential estimation procedure, as described in 
Appendix B, the state disturbance is viewed as white noise and is compensated 
for. An alternate method of handling state disturbance is to treat it as a signal 
considering its correlated character. Its expectation is still considered to be 
zero, but its autocovariance, considering the signal as a stochastic process, is 
associated with an assumed model. 

If the parameters of the auto covariance model are assumed to be known, 
then the signal and the state can be estimated as per normal least squares col- 
location procedures [Moritz, 1972], For example, consider the procedures in 
gravimetry. The parameters of normal gravity may be considered as consti- 
tuting the state . The gravity anomalies may be considered as the state disturb- 
ance or signal with known autocovariance function depending on the separation in 
distance instead of separation in time as in the case of this study. Further, the 
state would be considered time independent unlike the state in this study whose 
evolution in time is governed by the differential equations of motion. 

Butif the parameters of the autocovariance function arenotwell known, itis 
possible to include them inan augmented state vector in the adaptive filtering tech- 
nique or the dynamic model compensationtechnique, as described in [Myers, 1973], 

The detailed procedures depend on how it is proposed to model the state 
disturbance. Fo” example, in this study, referring to the statistics collected 
in Chapter 4, the state disturbance could be considered to be constituted mainly 
of position signal, velocity signal, or acceleration signal. Further, the 
acceleration signal could be either in the Cartesian coordinate system or in the 
spherical system. 

The computer software available was based on formulations for an 
acceleration signal in the Cartesian coordinate system. So, to judge the feasi- 
bility of using this procedure in the present problem, it was decided to attempt 
the procedure of adaptive filtering with the acceleration signal in the Cartesian 
coordinate system. The formulation for this based on [Myers, 1973] and 
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[Tapley, 1972] is given below. Formulations for other signals will need 
suitable modifications. 

The true equations of motion can be expressed as a system of first- 
order differential equations. 

r = v 

(5.2 

V = a (r, v, t) + i (t) 


where r and v are three-vectors which describe true position and velocity in an 
inertial frame as a function of time, a is a threc-vector function of acceleration 
components in the nominal dynamical model, and £ (t) is a three-vector accelera- 
tion signal. 

If £(t) is considered as white noise, the procedure of classical sequential 

estimation could be applied as stated earlier. But from the study in Chapter 4, 

it is indicated that the realizations of £-(t), £~(t), £^(t). which are 

the components of £(t), are correlated, and so a white noise model is not 

adequate. Therefore, £ (t) is approximated by a vector stochastic process 

7(1) (c (t) c (t) € (t)] T , and a simple model is assumed for this process 

X i z 

by considering the components of ?(t) as a time-correlated first-order Gauss- 
Markov process satisfying the stochastic differential equation 

t(t) - - B (t) 7 (t) + W € (t) (5.3) 

Here, W f (t) is a three vector of Gaussian noise with the properties 

E (W f (t) } 0 E i W f (t) WJ (s) j ---- q f (t) 6(t- s) (5. 4) 

where 

q f <t> 


o 

0 


0 

V 

0 


0 

0 




z 


(5. 5) 


is the matrix cf variance parameters associated with t (t), and B(t) is a 5 x 3 
diagonal matrix of time correlation coefficients 
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\ 0 

o 3v 


0 B r 


Further, for reasons explained overleaf, the evolution of 


Bit) = /J, 


is assumed to be governed by the differential equation 

S(t) - 0 ( 

Considering the above additional parametrization, the model in (5.2) is 
augmented to become 


r - v 

7 - a(r, v, t) + 7 (t) 

7 -1)7 + W (t) 

B o 


The new state vector is a 12-vector, 


X T (t) ( r T : V T : ? T : J§ T ] 


and the functional form of (5.8) is 


X = F(X, W, t); X (to) X 0 


(5. 10) 


With this augmentation the algorithm is similar to that of the first-order 
filter given in Appendix B. 

The advantage of the model of the form (considering one component). 


Lit) -a c {!)-■ \v f (t). 




(5.11) 


is that it gives rise to a stationary, exponentially correlated (laussian process 
or colored noise if the variance parameter and the time correlation coefficient 


are constant [Jazwinsky, 1970]. It may be recollected that an assumption of 
stationarity had been made for computing sample autocovariance of signals in 
Chapter 4. In this algorithm, the process is reinitialized at every step, and 
the parameters are treated as constant for that step. This explains why 
/3(t) = 0 in equation (5.7). 

From the above formulations, it is easy to understand the interrelations 
between some of the estimation procedures used in geodesy. In the adaptive 
filtering procedure used in this suidy, the state is time dependent. The state 
disturbance is viewed as a signal. The parameters of fee auto covariance 
model of the signal are assumed not to be perfectly known and are included 
in the sequential estimation procedure along with the signals. The state to 
be estimated would be given by equation (5. 9). 

If the parameters of fee autocovariance model of the signal are assumed 
to be perfectly known, these can be excluded from the estimation procedure, 
and fee state to be estimated would be 

X T (t) = [r T : v T : ? T ] 

This would be the familiar least squares collocation procedure extended to a 
dynamic situation. 

If the state disturbance is viewed only as a "noise " in this dynamic 
situation, fee state to be estimated would be 

X T (t) = [r T : v T ] 

as there is no signal to be computed, and this leads to fee Kalman filtering 
procedure or the first-order filtering procedure described in Appendix B. 

Further, if the state is assumed to be time independent in the above, 
the familiar least squares sequential adjustment procedure with weighted 
parameters is obtained. 
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5.3.2 Estimation of Initial Values of Parameters 


The sample autocovariances C££(u), C^,^.(u), C~ *(u) and sample 

autocorrelations r r - (u), r~-(u), r- "(u) are helpful to determine good initial 
XX YY Z Z 

values of parameters for use in the algorithm to improve the orbit represented 
by the broadcast ephemeris . 

The stochastic process considered above is the same as the first-order 
autoregressive process or the first-order Markov process [Jenkins, 1968]. 
Consider the component process c^(t) from above. It is assumed to satisfy 
a differential equation of the form 

'x ‘ -*x e x (t) + w 'x (,) ,5 ' 12) 


Considering the Gaussian white noise (t) as input and €^(t) as output, the 

X 

autocorrelation function of the output process, e (t), is 

X 


p xx (u) = 6 




[Jenkins and Watts, 1968, p. 162] 


(5. 13) 


where u is the time lag. This is an analog of the following discrete first-order 
autoregressive process 

= <514 > 


where 4* = E t e Y } . In this study 4=0, and so 
i A tt 


X., 


Ofi f Y + W 
X W t. 


(5.15) 


The autocorrelation function o f this discrete process is given by 

PxxW = a i k| (5. 16) 

where Qfi is the autoregression coefficient and k - 0, ± 1, ± 2 are the lags. In 
this study k ^ 0. The correspondence between (5. 13) and (5.16) is readily 
seen, which helps in computing from Ofi. 

Further, assuming a Gaussian distribution for W , the estimate ori of 

tfe 

ai is given by 
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[Jenkins and Watts, 1968] 


(5.17) 


A 




the autocorrelation at first lag (in this case, the lag of two minutes). 

It is clear that with taken from the r - -(u) curve in Chapter 4, 
the initial value of 0 parameter in the algorithm is easily computed. 

A 

Similarly, the variance parameter of W is estimated by 

ta 


s5 - c-wo, - 0l c--(i) 


for a large sample size. Parameters for other component processes and 
( are similarly obtained. 

Ci 


5. 3. 3 Simulation Program - EARTHOD 


The experiment in adaptive filtering has been carried out with the 
station positioning results of the experiment in Section 5.2, the statistical 
analysis shown in Chapter 4, and the simulation program, EARTHOD, 
obtained from the Department of Aerospace Engineering and Engineering 
Mechanics at the University of Texas at Austin, duly modified. 

The program provides the capability for simulation studies of an earth 
satellite observed by up to twelve ground-based tracking stations making as 
many as twelve simultaneous range, range-rate, elevation, and azimuth 
measurements. 

Observations are generated from a set of true equations of motion 
operated cm a true state of the satellite and are corrupted by Gaussian white 
noise. The program has an option for estimating the state of the satellite as 
well as the acceleration signal components and the parameters of autocovariance, 
based on the simulated observations, in an adaptive manner. 

Since both the true state and the estimated state based on the 
observations and the nominal state are available in the simulation, actual 
estimation errors are easily determined for the judgment of the filter per- 
formance. The equations of motion are expressed in an earth-centered 
inertial (EC!) coordinate system and are numerically integrated with an 
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efficient Runge-Kutta-Fehlberg algorithm [Fehlberg, 1968], 

For use in this study, the program written in Fortran IV in 
overlay form for the CDC 6600 operating system has been modified to run on 
the IBM 370/168 operating system at The Ohio State University, 

State vectors derived from the precise and broadcast epheme rides have 
been used in lieu of the initial true and nominal states, respectively, as inputs, 
after transforming them from the earth-fixed system to the inertial system. 
The propagation of state is carried out with an eighth-order Runge-Kutta 
algorithm. The station coordinates are assumed to be known. 

Results from EARTHOD simulations are available in terms of Root 


Sum Square (RSS) errors and covariance norms defined as follows: The RSS 
errors in position and velocity are, respectively: 


AR = l^ + 4 + e 2 z f 
AV = l4 + 4 +e Z^ 


(5. 18) 

(5. 19) 


where e ~ X - X represents the true error component in the state position 

X A 

estimate X. Similar definitions apply to other components, e • = X - X 

A A 

represents the true error component in the state velocity estimate X . 

In this study, X represents the vector of state components based on 

A 

precise ephemeris and 3T represents the vector of their estimates based on 
observations and the nominal state component. The nominal state was 
derived from the broadcast ephemeris. 

The position covariance norm is defined 

i 

NP - tPn + Psb + Psb] (5.20) 


the square root of trace of the covariance elements associated with the posi- 
tion estimate and is obtained from the diagonal of f* matrix in the algorithm, 
vide equation (BJ>8). 

The velocity covariance norm is dofined by 

k 

NV fp* + Pi*+Pi«r (5.21) 
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which is likewise the square root trace of the covariance elements associated 
with the velocity estimate. 

5. 3. 4 Description of the Experiment and Results 

An experiment in adaptive filtering was carried out for one 
pass of satellite 19 in Data Set D. The statical coordinates for the 
three stations, as obtained in the experiment (before removal of radial bias), 
were taken as known. Within the capability of the software, range-rate obser- 
vations were the closest to Doppler observations, and so for further reasons 
explained in Section 5.1 the mode of simulated range rates was chosen, giving 
a three-vector observation at each step. 

A fixed integration step of 0.075 minutes was chosen which is close to 
the rate of data acquisition with a JMR 1 receiver. The common two-minute 
epoch at the commencement of the observations from the three stations was 
chosen as the initial time t 0 for commencing the algorithm. The precise state 
vector (transformed to WGS 72 system) and broadcast ephemeris state vector 
at to were transformed from the earth-fixed system to the inertial system. 

The values for polar motion components were taken as published by Bureau 
International de l'Heure, and the value for the Greenwich apparent sidereal 
time at to was obtained from the American Ephemeris and Nautical Almanac 
1976. The Jet Propulsion Laboratory ephemeris tape provided the inertial 
state vectors for sun and moon. 

In the experiment for obtaining station positions, the range observation 
standard deviation was taken as 0.20 m. So a range- rate standard deviation of 
0.002 m/s appears reasonable. Based on the study in Chapter 3, the uncer- 
tainties in the position and velocity components of the broadcast state vector 
were taken as [ 14.0 15.0 22.0 0.20 0.25 0.30 ] in units of meters and 

A 

seconds for initial values in Po in the first six diagonal locations. The values 

for the next three diagonal elements in Po, which referred to the autocovariancc 

of the acceleration signal at lag 0, were taken from the curves C- - (u), C- - (u), 

AA Y Y 
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and C~ -(u). The last three elements of P 0 refer to the uncertainties of the 

4 z 

time correlation coefficients. The values for these were taken either as five 

times the coefficients (as used in past investigations of Myers [1973]) or 

derived from the difference between the two values for the coefficients from 

the two sets of the curves C-~(u), C - " (u), and C~ - (u) for satellite 19, one 

XX YY Z Z 

from Data Set D and the other from Data Set S. 

Four trials were taken, varying slightly the standard deviation of 
observation and the initial values for the state noise parameters, from the 
analysis in Chapter 4. 

In trial 1, the standard deviation of observations was kept as 0. 002 m/s, 
but the initial parameter values were based on the information collected for 

A 

satellite 19 from Data Set D. Uncertainties of fi , fi. fi„ for P 0 matrix 

x y ^ 

were taken as five times their values. 

In trial 2, the standard deviation of observations was kept as 0. 002 

m/s for the R matrix In the algorithm, but the initial parameter values for 

state noise referred to Data Set S. Uncertainties of fi , 8, and fi for 

x y z 

A 

P 0 matrix were taken as five times the initial values of the coefficients . 

In trial 3, the standard deviation of observations was increased to 0.005 
m/s, and the initial parameters were the same as in trial 2 . 

Trial 4 was similar to trial 2 in all respects except that the uncertainties 

A 

of fit fi , fi in the P 0 matrix were based on the difference between the two 
x y & 

sample values for fi , fi. fi for satellite 19 obtainable from Data Sets D 
x y z 

and S. 

The structures of the A and H matrices required in the algorithm are 
given in detail in [Ingram, 1971] and [Myers, 1973] and incorporated in the 
software used. 

The sample autocovariance curves from the data analyzed in Chapter 4 
were used directly for the trial values with the procedures given above. No 
least squares adjustments were made for estimates of parameters since it 
was the intention of this experiment to judge the adaptive property of the filter. 


89 


Figs. 5.1 and 5. 2 show die values of AR, NR and AV, NY, respec- 
tively, against i, the number of integration steps of 0.075 minute, for trial 
1. Figs. 5. 3 to 5. 8 give the same information for trials 2, 3, and 4. 

Table 5-3 gives the components of state vector pertaining to the state 
disturbance (i.e., acceleration signal and time correlation coefficients) and 
their values at various integration steps, i. Tables 5-4 to 5-6 give the same 
information for trials 2, 3, and 4 . 

5.3.5 Observations 

From the results of die trials, the following observations can be made: 

(i) The initial deviation of 12.91 m in die position of the satellite as given by 
precise and broadcast ephemerides narrows down by about 0 . 6 m in 
about 20 integration steps but later diverges . 

(ii) The uncertainty in position of the satellite as given by broadcast ephemeris 
reduces by about 4. 2 m in about ten integration steps and later diverges. 

(iii) The rate of change of improvement is the largest during the first five 
integration steps. 

(iv) The optimum situation for position uncertainty arises fairly close to the 
optimum situation for position and die AR curve is flatter than the NR 
curve. This is of particular significance from the point of view of users 
who have no access to precise ephemeris. When this procedure is used 
with real observations and the broadcast ephemeris, only the NR curves 
are available . The curves in this study indicate that the state vector at 
optimum positional uncertainty, given by the NR curve, also yields a near 
optimum value for the position of the state vector. 

(v) Velocity uncertainty decreases by about 0.2 m/s in the first five integra- 
tion steps. This is followed by further improvement at a very slow rate. 

(vi) Velocity error reduces by about 0, 01 1 m/s in the first step. Later the 
error increases steadily to about 0.033 m/s in about thirty integration 
steps. This is followed by further steady improvement. 
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Figure 5.1 

Position Error and Covariance Norm in Trial 1 
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Figure 5.2 

Velocity Error and Covariance Norm in Trial 1 










Figure 5.6 
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Table 5-5 


State Disturbance Components in Trial 3 


Integration 

~ri 


n 

1 

sec 

X 

sec 

l /* 

IK 

0 

0. 00 

0.00 

0. 00 

163.20 

wm 

1 

5 

-0.0 1 

0. 10 

0.02 

168.20 


■ 

10 

0.00 

0.03 

-0. 00 

163. 20 

81.48 


15 

0.02 

-0.03 

-0.07 

168.2 0 

8 1.4 0 

134.58 

20 

0.04 

-0.07 

-0. 13 

168 . 22 

31.05 

134.70 

25 

0.06 

-0.11 

-0.20 

168.31 

82. IS 

135 . 18 

30 

0.09 

-0. 13 

-0.27 

168.49 

82.32 

136.22 


0 . 10 

-0 . 14 

-0.34 

168 . 80 

83.34 

137. 95 


0. 12 

“0. 13 

I 

-0.40 

169 . 20 

8 3-60 

140. 28 

wm 

0, 13 

-0.11 

-0. 44 

163.63 

8 3. 52 

14 2 . 97 

50 

0. 13 

-o.os 

-0. 46 

1 70. 04 

83. 16 

145. 67 


Table 5-6 

State Disturbance Components in Trial 4 


Integration 

step 

i 



■ 

4 

sec 

H 

sec 

0 

l) . 00 

0 . 00 

0.00 

168 . 20 

3 1.60 

134 . 59 

5 

0 . 00 

-0.03 

0 . 02 

166. 20 

8 1.60 

134 59 

10 

0.01 

-0.16 

-0.02 

168 . 20 

8 1.60 

t34.59 

15 


-0.26 

-0 .02 

168.20 

8 ! . 60 

134.59 

20 


-0. 34 

-0.07 

168 . 20 

8 1.60 

134.59 

25 

0.1) 

-0.38 

-0.13 

168 . 20 

i 1 . 6 1 

134.59 

30 

0 . is 

-0.41 

•0 , 19 

163. 2 J 

8 1.62 

134.59 

35 

0.18 

-0.40 


168 . 22 

8 1.62 

134 . 59 

4 ) 

0 . 20 

-0.36 

-0.31 

168 . 22 

3 1,62 

134 . 60 

45 

E ■ 

-0. 2 9 

-0.36 

168 . 23 

3 1 . «i 2 

134.60 

50 


-0.22 

- 0 .3 9 

168. 24 

8 1.62 
i— 

134.61 
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(vii) From the observations at (i) to (vi), it spears that from the point of 
view of a user, the value of the state vector at optimum uncertainty of 
position could be taken as the best outcome of this filtering procedure. 
Specifically, in this experiment the results of the tenth integration step 
would be the optimum . 

(viii) From Tables 5-3 to 5-6 it can be seat that the values of 0 X> ( 3 ^, 0 Z 
are relatively insensitive to the filtering procedure. So it is best to 
use as realistic values as possible for this procedure . On the other 
hand, changes in values of 0 » 0 , 0 in trials 1 and 2 have not given 
significantly different results. 

This indicates that while it is advisable to repeat studies of this nature 
occasionally to ascertain the autocovariances for all satellites, since 
the autocovariance curves of the three satellites in this study are fairly 
close to each other, the values based on the data analysis in Chapter 4 
could be applied to other satellites as well, until studies for 
other satellites are completed . 

(ix) It is clear that some improvement in position and velocity is possible 
with this procedure though the gain in reducing the uncertainty of position 
and velocity is greater. It is, therefore, conceivable that if the software 
is suitably modified the optimum state vectors for the 12 passes in experi- 
ment 1 can be obtained using broadcast ephemeris and observational data. 
With these optimized state vectors held fixed, the solution in experiment 1 
could be repeated for improved station position recovery. 

(x) It is noticed that the filter diverges within a very short time of about 20 
integration steps. However, from Chapter 4 it is known that this filtering 
procedure is being used with acceleration signal in Cartesian coordinate 
system from consideration of available software. This signal has a posi- 
tive autocovariance for only two to four minutes, and the signal compo- 
nents are mutually correlated to some extent. Much better results can 
be expected if the positional signal in spherical system is used, which 
has a positive autocovariancc up to about 20 minutes. However, the 
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feasibility of this procedure has been demonstrated. The divergence is 
also partly due to the weak geometry represented by three nearby 
stations (onlyabout30 km apart), coobserving a satellite at a height of 
about 1100 km. 

(xi) As expected, the results of trial 3 with a larger observational standard 

deviation show a smaller improvement compared to the other trials. 

(xii) Results of trial 4 are the best among the four trials. 

(xiii) The signals e , c , f can be viewed as acceleration biases associ- 
a Y Z 

ated with the specific pass of the satellite used in the experiment. It 
can be easily seen that if a filtering procedure as above is carried out 
with formulations for positional signals in a spherical system, the 
signals obtained would be the in-track, out -of -plane, and radial pass 
biases referred to in Section 3.5.2. 

5.3.6 Limitations of the Filter World 

For computing A R, the positional deviation between the estimated state 
and the true state of the satellite, the latter is obtained by integrating the state 
derived from the precise ephemeris at to. Ideally, the force model for this 
integration should be the one adopted by DMA with all its elaborations. Simi- 
larly, the nominal state should use the force model of the NAG computations. 

However, due to the limitations of the software and the nonavailability' 
of the values of the gravity field in current use by DMA and NAG in open 
literature, there is an inescapable mismatch between the force model in the 
filter world and reality. 

The force model in the filter in this experiment includes the gravitational 
field of the earth (GEM 7 geopotential model [Lerch et al . , 1975) with coefficients 
up to degree and order four and additional zonals up to degree six), two-body 
perturbations of the sun and the moon, and atmospheric drag, for both the 
true and the nominal state. 
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So in order to obtain realistic results , it has to be assumed that the 
optimum results of the filter are obtained before they are vitiated by the 
mismatch in the force models . 

For checking the validity of this assumption, trial 4 was repeated 
with an integration interval of 0. 05 minutes, and the positional deviation of the 
estimate was computed directly with respect to the corresponding precise 
ephemeris state vectors at integral minutes. Fig. 5.9 shows this deviation 
AR P against the time t in intervals of one minute after to. 



Figure 5.9 

Position Error with respect to Precise Ephemeris in Trial 4 


Comparing this curve with Fig. 5.7, it can be seen that there has been 
no adverse effect of the mismatch during the first minute in which the optimum 
results of the filter have been realized. The three-minute mark in Fig. 5.9 
corresponds to die 40th integration step in Fig. 5. 7. 
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6. CONCLUSIONS AND RECOMMENDATIONS 


The objective of this study has been to investigate the possibility of 
locally optimizing the Navy Navigation Satellite broadcast ephemeris for 
improved recovery of station positions with Doppler observations. 

The rank deficiency problem which is intimately connected with Doppler 
surveys has been studied, and it has been shown that the minimum rank 
deficiency in a short arc mode of survey is six and the scale information is 
derived from the wavelength of transmission. Coordinate differences are 
estimable quantities if the stations coobserve the same pass of the satellite 
whose motion is governed by an assumed force model. 

Accuracy estimates of broadcast ephemeris have been formed from the 
study of sampled data . It is concluded that, depending on the location of the 
epoch of obse iwat ion in the interinjection period, the positional uncertainty of 
broadcast ephemeris may vary between 19 m to 26 m in-track, 15 m 
to 20 m cross-ti-ack, and 9 m to 10 m in radial directions. 

The broadcast ephemeris indicated a radial bias of about -5 m when 
compared with the precise ephemeris transformed to WGS 72 system accord- 
ing to the currently known transformation parameters . 

The experiment in removal of radial bias along with recovery of station 
positions holding the orbit fixed showed that the removal of bias has a negli- 
gible effect on the uncertainties of the station coordinates. But the values of 
the coordinates as well as the values of other pass parameters change slightly. 
Removal of bias brought a chord distance in better agreement with terrestrial 
survey. This radial bias appears to be the consequence of the following two 
factors brought to light in this study; 

(i) The origin of the broadcast ephemeris coordinate system appears to 
have an offset of about 10 m in the direction of 90°W longitude, with 
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respect to the geocenter. 

(ii) The broadcast ephemeris appears to need a scale correction of -1. 4 ppm 
to make it compatible with a terrestrial system obtainable from a scale- 
corrected precise ephemeris. 

It is recommended that experiments be carried out with more data sets 
in the future to identify in-track and out-of-plane biases also. These could not 
be included in this study for lack of adequate evidence with the available data 
sets. 

Sample autocovariance and autocorrelation functions of acceleration, 
velocity, and position signals have been computed by comparing precise and 
broadcast ephemerides. The curves indicate an exponential form which has 
been assumed for the signal stochastic process. 

Satellites decay with time . So it is recommended that studies similar 
to this study be repeated occasionally so that the statistical characteristics 
of the broadcast ephemeris of all the satellites are available as they evolve 
with time. It is difficult to suggest the time interval at which such studies are 
to be repeated. In this study two data sets of the same satellite (satellite 19) 
were available at an interval of ten months . But the earlier data set had only 
about 180 values against about 1500 in the second. The same autocovariance 
functions derived from the two data sets do show variations which are partly 
due to the unequal size of the data sets . Keeping this in view, an interval of 
one year for repeating such studies appears reasonable. 

The experiment in adaptive filtering for one pass of satellite 19, with 
station coordinates obtained in the earlier experiment, parameters of the auto- 
covariance model of acceleration signal in the Cartesian coordinate system 
obtained from sampled data, and simulated range rate observations, leads to 
the conclusion that adaptive filtering is a feasible approach from a practical 
point of view. The experiment indicated an orbit improvement of about 4 m in 
the positional uncertainty and 0.6 in position within the first 20 sequential co- 
observations from three-stations at intervals of 0.75 minutes. But the filter 
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diverges after giving this improvement. Positional signals in the polar system 
which have autocovariances which are larger and remain positive for up to 
about 20-minute lag are expected to give better results with this procedure. It 
is recommended that further studies in this respect be made. 

Any satellite ephemeris derived from satellite transmission is of 
necessity predictive and less accurate than that of an ephemeris computed 
after the fact and possibly with more data (e.g. , the satellites of the NAVSTAR 
Global Positioning System) . In all such situations where sample statistics of 
state disturbance can be predetermined from comparison of predicted and 
post-fitted ephemeris, filtering procedures could be applied in conjunction 
with broadcast ephemeris and current observational data for local improve- 
ment of orbit. 

Further studies are also recommended to study the performance of 
the adaptive filter with real data and a better geometry of tracking stations 
than three stations about 30 km apart as in this study. 

As a result of this study, the following steps can be recommended for 
improved positioning with broadcast ephemeris and Doppler data: 

(i) Correction of broadcast ephemeris for the biases identified. In this study 
the radial ' .as identified was traced to an origin shift and a scale correc- 
tion. 

(ii) Position determination with broadcast ephemeris and coobserved tracking 
data from at least three stations, holding the orbit fixed. 

(iii) Local improvement of the satellite orbits for all the passes used in (ii), in 
an adaptive filtering procedure, using the sample autocovariance function 
derived from comparisons of precise and broadcast ephemerides in the 
recent past and the observational data. 

(iv) Redetermination of positions as in (ii) with the improved orbits. 

The improvement that can be achieved by this procedure can only be 
estimated after studying a test case with real data for steps (iii) and (iv) which 
is recommended as a follow-up of this study. However, this study has 
succeeded in indicating the feasibility of this approach. 


102 



APPENDIX A 


INFLUENCE OF TRANSFORMATION PARAMETERS 
ON ACCELERATION SIGNAL 

The broadcast ephemeris and the precise ephemeris are in coordinate 
systems which differ from each other by small amounts. The following ana- 
lytical study was carried out to find out to what extent the transformation 
parameters are likely to influence the values of acceleration signals obtained 
in Chapter 4 . 

Consider the following system of first-order differential equations for 
governing the true state of the satellite in the WGS 72 system of the broadcast 
ephemeris. 

= Vn (Al) 

Vn = a N + 4 n (A 2) 

where X N , V N , a N , £ N are the three-vectors for position, velocity, acceleration, 
and acceleration signal, respectively. 

Similarly, consider the following system of equations for governing the 
state of the satellite in the NWL 9D system of the precise ephemeris, which for 
the purpose of this study is considered errorless: 

= V p (A3) 

V P * a p (A 4) 

where X P , Tp, a P denote the three-vectors for position, velocity and acceleration 
as per precise ephemeris. 

It is assumed that the state vectors based on the precise and broadcast 
ephemeris represent a set of consistent values in their respective coordinate 
systems, which makes the following similarity transformation relation possible 
for the i^ 1 point 
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(A 5) 


Here 

[X p Yp Z/ = X p 

[X N Y n Zn] t = X N 

[ AXp AYp AZ P ] T = T p , the vector of translation parameters . 

A L is the scale parameter 

W, e, and $ are the rotation parameters. 


If R is the 3x3 rotation cum scale matrix, equation (A 5) can be rewritten as 


X N = T p + RX P (A 6) 

The transformation parameters are time independent, so time derivatives of 
(A 6) yield 

X N = RXp+RXp = RX P as R = 0 (A 7) 

Therefore, 

V N = RVp (A 8) 

V N = RV P + RV P = RV P (A 9) 

a N +I N = R[a p ] (A 10) 


To seek the influence of transformation parameters on £ * signal, consider 
a vector 

U = [AL, W, *, € ] T 

and differentiate (A 10) with respect to the variables of interest. The left 
member yields, on differentiation, 
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considering R matrix in the element notation, the right member of (A 10) may 
be written 


R[a p ] 


" r n 

r 12 

ria" 


r n a P i 


ria a w ‘ 


ri3 aja 

r a 

rae 

r» 

„ P| 

II 

r a a rt 

+ 

rae && 

+ 

r SB a* 

_r 3 i 

raa 

T33_ 


r 3 i a pi 


ras aps 


1*33 a ]8 


Ri + R 2 + Rs 


where R lf R 2 , Ra stand for the corresponding vectors above. Differentiation 
of the right member therefore yields 


If dC» + ^ du + ^5dC, + ^ dy + ^ dC» + — du 
3U SC, 81? »?, au 


The modeled acceleration a N does not depend on the acceleration signal or on 
the transformation parameters . Therefore 



Similarly, 

d Ri _ d Rg _ dRa 


Therefore, differentiation of right member of (A 10) yields 


fa* du + ff du + |f du = 
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Therefore, 


d$, 


j .111 . 4 MW MI M W^ Wt iQWtC y-.^f^ 41 


9 rn 

9AL 

9r n 
9 \V 

dr a 

9AL 

9 r a 
9 VV 

9rai 

9AL 

i£a 

9 W 

9 ria 

9AL 

9 ri? 
9 W 

9 raa 

9AL 

9rg» 
9 W 

9 r^ 
9AL 

9 rag 
9 \V 

9 r u 
9AL 

9 ru 
9VV 

9 ra 
9AL 

ila 

9\V 

9 r^B 

9AL 

9 ras 
9W 


iin 

9* 


3 Ti 
9c 


ila 

9 


9 SI 


9c 


ila 

9* 


9 r ; 
9c 


flpi dU + 


3 r 
9 ^ 


9ria 
9 c 


9 ra: 
9* 


9 r a 
9c 
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9* 


9 r ; 
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5 r^ 
9 ^ 


9 r 
9c 


9 r g 

9* 


9r^ 
9 c 


9 rg; 
9 ^ 


9 ras 
9c 


3ps dU 


1 

0 

0 

1 

o 


0 

1 

0 

o 

0 

-1 

0 

0 

a P idU + 

1 

0 

0 

0 

0 

0 

1 

0 


0 

0 

0 

-1 


a P g dU + 


0 -1 

0 0 


0 


a p3 dU 


0 

oj 



ft pi 

ftp? 

~ftp3 

0 

ft pa 

~ftpl 

0 

a p3 

ft p3 

0 

ftpl 

~ a pi? 


dU 


(All) 


In this study, the U vector has the following values for the components 
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AL = -0.826 3 ppm 

W = -0.000 001 260 5 radian 

$ = 0 

e =o 

At the satellite altitude of about 1000 km, the approximate value of 
normal gravity is 743.4 gals. Considering this value and even the full values of 
the U vector, the maximum component of d£ N is of the order 0. 15 x 10“ 6 m/s 2 , 
which is negligible compared to the signals of the order of 0.3 x 10” 3 m/s 2 
obtained in the study. 

It may, therefore, be concluded that the transformation parameters have 
had an insignificant role in causing the acceleration signals. 

This was also verified numerically by obtaining acceleration signals in 
two models in a test case: 

(i) by comparing broadcast and precise ephemeris without transformation, 

(ii) by comparing broadcast and precise ephemeris after transformation. 
Considering acceleration in units of m/s 2 , the RMS values of the signals in (i) 
and (ii) are generally in agreement up to the fifth decimal place of a meter. 
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APPENDIX B 


FIRST ORDER FILTERING TECHNIQUE 


1. Systems Equation 

The true equations of motion for the satellite are expressed in the earth- 
centered inertial coordinate system. The equations take the functional form of 
a system of first-order deterministic differential equations (i.e. t without the 
inclusion of state noise): 

r - v 

(Bl) 

V ^ a(r, v, t) 

where Twith components (X, Y, Z) and v with components (X, Y, Z) are three 
vectors used in defining the true state X, an n-vector which is to be estimated. 
For example, if the state vector were to consist of only position and velocity 
components, then n - 6, and 

3T T - [X Y Z X Y Z) T = [Y t V t ] 

a is the true acceleration vector. Equation (Bl) can also be written in a 
functional form as 

* F(X, t) (B 2) 

where F represents the n-dimensional functional form of the right side of 
the equations (B 1). 
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Now, if a state noise n-vector W (t) is also included in the equation, the 
differential equation will take the form 

t F ( X , W, t) (B 3) 

For seeking a solution, a nominal state X *(t) is first assumed, and a differ- 
ential state deviation is considered in the form: 

x(t) •-= X (t) - X*(t) (B4) 

j(t) - ?t(t) - %*(t) (B 5) 

State relation (B3) can now be expanded by a Taylor series around X*(t) giving 
the linearized first-order differential equation 

j(t) = A (t) x(t) + G(t)W(t) (B 6) 

which corresponds to the equation given in (4.1) [Ingram, 1971; Tapley et al. , 
1972]. Here 



are the n x n matrices of partial derivatives evaluated at the nominal values 
X*(t). 

Equation (B6) represents the state of the system (in this case the param- 
eters associated with the satellite orbit) for t > to, where the initial time t 0 is 
fixed and the initial state X (to) is assumed known. 

2 . Concept of the State Transition Matrix 

In many problems a representation in terms of the so-called state 
transition matrix helps to obtain an explicit expression for the solution of 
equation (B6). Since the system in (B6) is linear, its complete solution 
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consists of a linear combination of its homogeneous (or force-free solution) 
and a particular (or forced) solution. 

Consider first the homogeneous solution of the system 

j(t) = A(t) x (t) (B 7) 

for t ^ to with x (to) arbitrary. Substitute into (B 7) a trial solution of the 
form 

x(t) = M(t) x(tb) 

where M(t) is an unknown n x n matrix. This yields 

[ M - A(t) M ] x (to) = 0 (B8) 

which must hold for all t s to. Since 3? (to) is arbitrary, this relation is satis- 
fied if and only if M(t) satisfies the system of n x n differential equations 

M = A(t) M (B9) 

for all t ^ to. Ptirther, at t = to, x(to) = M(to) x(to) which implies that 

M(to) = I is the initial condition for (B9). The homogeneous solution of equation 

(B6) can then be written as 

*(t) = M(t)x(to) (BIO) 

where 

M = A(t) M and M(fo) = I 

The particular solution for (B 6) is next obtained by using the Lagrange variation 
of parameter technique. Assume a solution of the form 

Jf (t) = M(t)T(t) (B 11) 

where M(t) is as above, and 1 (t) is an unknown n vector. Substituting this result 
in (B6) 

M(t)7(t) + M(t)7(t) = A(t) M(t)7(t) + G(t)W(t) (B 12) 

However, since M(t) = A(t) M(t), (B 12) reduces to 

M(t)7(t) - G(t)W(t) (B 13) 
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It can be verified that M(t) is nonsingular for all t s t 0 [IV.'editeh, 1969] . 


Therefore, 


1 (t) = M x (t) G(t) W (t) 


(B 14) 


Integrating (B14), 


1 (t) =1 M 1 (T) G(r) W(r) dr 

tc 


(B 15) 


Therefore, the particular solution of (B 6) can be taken as 


x(t) = M(t) J M 1 (T) G(r) VV(7) dT 


(B 16) 


Combining (BIO) and (B16), the complete solution of (B6) will have the form 


x(t) = M(t) x (to) + M(t) C M 1 (T) G(T) W(T) dT 

tb 


(B 17) 


where the second part of the right member can be viewed as the contribution 
of the state noise to the state. 


M(t) is called the fundamental matrix of the system in (B6). Now 


define 


$(t, T) = M(t) M 1 (T ) 


(B 18) 


the n x n transition matrix for the system in (BG). It is noted that 
*(t, to) - M(t) M 1 (to) - M(t) 


since M(to) - I. Thus the complete solution (B 17) can be written in transition 


matrix form as 


x(t) -- $(t, t 0 ) x(to) + J $(t, T)G(T) W(T) dT 


(B 19) 


where t - to. (B 18 ) is differentiated with respect to t 
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= M(t) M l (T) = A(t) M(t) M'V ) 

Therefore, 

$(t, T) = A(L)*(l t T) (B 20) 

where differentiation with respect to t is implied. From (B 18) 

$(t, t ) = I, for t ^ to. 

So, to obtain $ (t, T) required in the computation of x (t) in (B 19), the 
differential equation 

4(t, t 0 ) = A(t)$(t, to) 

is solved subject to initial condition $(to, to) = I and to is replaced by T to 
obtain ^(t, T). 


3. Qbseryp ticn - State Relation 

The equation relating observations and the state in a linearized obser- 
vation equation form can be expressed as 

y(t) = H(t) 5E(t) + V (t) (B21) 

where y(t) is the p vector of observations in the differential form (observed - 
computed). H is the p X n design matrix of the conventional observation 
equation model [Uotila, 1967], x (t), the deviation of thestate from the nominal, 
can also be viewed as the unknown corrections to the assumed values of the 
state. V (t) is the p vector of observational errors . 

Substituting the expression for x(t) from (B 19), 

t 

y (t) = H(t) [$(t, to) x(to )+f $(t, T) G(T) W(T) dr] + V(t) 

to 

t 

H(t) <t»(t, to)x(to) + J H(t)*(t, T) G(T) W(T) dT + V(t) 
to 
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Defining 


D(t, T) = H(t)$(t, T) G(T), 

y (t) = H(t) $(t, to)3T(to) + j D(t, T) W(T) dr + V(t) (B 22) 

to 

Often the noise vector W(t) is considered as input, the observable y (t) as 
output, and the matrix D(t, T) as a system weighting function, in the lit- 
erature. 


4. Sampled Data Model 


Although the evolution of state is continuous in time, observations are 
generally available only at discrete times. The study of the effect of noise is 
also facilitated if the system is discretized and the continuous relations in a 
system are treated as limiting cases of discrete forms. 

Assuming that the disturbance or "noise" vector in equation (B6) is 
a piece-w'ise constant function of time which changes values at time points at 
which measurements are also made, a time interval t k £ t s t k +i is considered 
for some k ~ 0, 1, ... . If x(t k ) is given and W(t) = W(k) - constant, in the 
interval t k s t s t k +i, from (B 19): 


X(tk + l) 


$(tk+l, t k ) x(t k ) 



r) G(r) dr 


W(k) 


(B23) 


Defining 

tk+l 

f ^ (tk+i» r ) G(r ) dT - r {k+ 1, k) 

tic 

as the disturbance transition matrix, and denoting 

x(tk+i) x(k + 1) 

x(t k ) 5c (k) 

*(tkn, tk) *(k+l, k) 
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(B 23) can be written as 

x(k+ 1) - $(k + l, k) x(k) + T(k + 1, k) W(k) (B24) 

for k ' 0, 1, ... . 

The observation state relation (B21) can also be written in the 
discretized form as 

y(tk+i) ^ H(t* +1 )x(tk + i) + V(t k+1 ) 
or 

y (k+ 1) * H(k+1) x(k+l) + V(k + 1) (B25) 

It is easy to see that the continuous case can be considered as a limiting case 
of the above discrete model by denoting discrete time (k) and (k+ 1) as (t) 
and (t + At) and letting At -* 0. 

5 . Nature of Gaussian White Process 

The nature of Gaussian white process in system dynamics will be 
first considered. Let {\V(t), t s to } be an n-dimensional independent 
Gaussian process with mean 

E [ W(t) ] W„(t) 

and covariance kernel 

E l [W(t) - W»(t)] [ W(T) - W B (T )] T } - Q(t) 6(t - T) 

where to is an initial time, t, T ^ to, Q(t) is a continuous positive semi- 
definite n x n matrix, and 6(t - T) the Dirac delta function. 

The above equation implies the limiting behavior of a piece-wise 
constant Gaussian white sequence in which the frequency of event points is 
made arbitrarily large within a given time interval. 

Consider sequences fW(k), k 0, 1, ... ] to be zero mean Gaussian 
white with covariance 

E [ W(j) W T (k)l Q(k) 6 Jk , k 0, 1, ... 


114 



and successive time points separated by At > 0. If t denotes continuous time, 
to corresponds to k = 0, and ti corresponds to k = n, 

ti = to + n A t 

For given value of n, let { W (n) (t), to s t ti } denote a piece-wise constant 
Gaussian white sequence. Keeping ti fixed while increasing n such that 
nAt = ti - to is constant, the nature of { W Ca) (t), to 5 t ^ ti } as n 00 and 
At •* 0 is considered. 

As proved in [Meditch, 1969], the Gaussian white process 

(W(t), to* t* t! } = lim tW^ft), 1o s ts ti } 

n~> 00 

as described above, where the covariance matrix Q(to + i At) = Q(i), i = 0, 1, 

. . . , n-1, is to be replaced by Q(t)/At in taking limits where t corresponds to 
the time point i and Q(t) = Q(i) . Then ^ is understood in the sense 

that the quantity dealt with is defined over an interval of width At, and in the 
limit as At - * 0, the function which isl/At over the interval At and zero 
elsewhere becomes the dirac delta function, giving 

E C [W(t) - W.(t)] [ W(T) - W b (T)] T } = Q(t) 6(t - T) 

for all t, T 2 to. 

6 . Probabilistic Description of System Dynamics 

The probabilistic nature of state noise and its constribution to the 
evolution of the state makes the study of probabilistic description of the state 
model imperative. 

6.1 Evolution of State in Presence of White Noise 

As indicated earlier, equation (1124) can be rewritten in the form 
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t+At 


x(t- 


At) = $(t+At, t) x(t) + f $(t+At, T) G(T) W(T) dT = 


= 4>(t+At, t) x(t) + T(t+At, t) W(t) 


(B 26) 


under the assumption that W(T) - \V(t) - constant for t 5 T < t + At. Treating 
the noise process I W (T ), t s T ^ t + At] to be the limit of a sequence, as 
described above, W(t) in (B26) can be replaced by a member of the sequence 
W Cn) (t). Further, ^(t+At, t) expanded around t in a Taylor series up to 
linear terms yields 

$(t + At, t) $(t, t) + i(t, t) At = I + A(t) At (B 27) 


Als o. 


r(t + At, t) 


t+At 



<t>(t+At, T) G(T) d T 


G(t) At 


(B28) 


With these simplifications, equation (B26) can be written as 

x(t + At) = [I + A(t) At] x(t) + G(t) W <n) (t) At 

This can further be written in die form of a difference equation as 

x(t + At) - x(t) A(t)x(t) At + G(t) W (n) (t) At (B 29) 

Now letW(t) be a Gaussian white process with the mean E { W(t) } - W B (t) 
and covariance kernel 


E (f W(t) - W«(t)] [ W(T) - W,(T )] T } Q(t) 6(t - T) (B 30) 

where Q(t) is a continuous semidefinitc n x n matrix, and 6(t - T) the dirac 
delta function. Now dividing (B29) by At, and taking the limit At > 0, leads 
to the continuous linear system 

£ - A(t) x •* G(t) W(t) (B 31) 

Now let x (to) be a Gaussian random n vector which is independent of { W(t), 
t ^ to ) nnd has mean x B (to) and positive semidefinite n x n covariance matrix 
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E { [ x (to) - x B (to)] [ x (to) - x.(to)] T } = P(to) 

As a consequence of the assumption of independence, 

E { [ x(t 0 ) - Xc(to)] [ W(t) - W m (t)] T } = 0 

for all t ^ to. 


(B32) 


(B33) 


To examine the nature of ( x(t), t ^ to }, the solution of (B31) can be written as 


x(t„) = $(t», t rl )x(t rl ) + j *(t„ T)G(T)W(T)dT (B34) 

trn- i 


Thus, for t 0 > tu-x - to, the probability law describing the process x(t) in the 
future (i.e., at time t„) depends only on the present value the process assumes 
(i.e., at tn-i) and is completely independent of the behavior of the process in 
the past. Therefore, the process x(t) is a Markov process. 

Further, if x(to) is a Gaussian random n vector, then it can be 
deduced that 

x (a) (t), n = l, 2, ... 

is also a Gaussian random vector [Med itch, 1969], Thus x(t) is a Gauss- 
Markov process. Taking expectations of (B31), 

X, = E {£} A(t) x B + G(t) W B (t) (B 35) 

for t ^ to and subject to initial condition x (to) . 


6 .2 Evolution of State Covariance Matrix in Presence of White Noise 


An expression is now obtained for 
P(t) E ((X(t) - x B (t)] [x(t) - x.(t)] T } 

From (B 26) 

x (k + 1 ) 4> (k + 1 , k) x (k) + T (k ^ 1 , k) W(k) (B23) 

from which it is clear that 

x(i < 1) 4>(i 4 1, i) x(j) 4 r (j * 1, j) W(i) (B 36) 
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and 

x(j + 2) = *(j + 2, j+1) x(j + l) + TO + 2, j + 1) W(j+1) 


Therefore, 

x(j+2) 


= ' *0 + 2, j+1) [ *0 +1, j) x(j) + ro + 1. j) W(j) ] + 

+ r (j + 2, j + i) wo+i) 




*0 + 2, j) xfl) + ^ $ (3 + 2. i) r (i, i-1) W(i-l) 


J+1 


(B 37) 


Following the above recursive x’elation, 

k 

x(k) = * (k, j)x(j) + *(k, i) T(i, i-l)W(i-l) (B 38) 

i= j+i 

It is clear that for W(i - 1), i ~ 1, . . k, Gaussian, x(k) is Gauss-Markov. 

From (B26), taking expectations, 

x„(k+ 1) --- $>(k + l, k)x B (k) + T(k+ 1), k) W B (k) (B39) 

for k € I, the set of integers. If x B (0) and W B (k), k € i are given, (B 39) 
becomes a recursive relation for the mean of the sequence. Now 

P(k+ 1) - E {[x(k+l) - x«(k + l)nx(k + l) - x,(k + l)) T } 

K ( t*(k+l, k) [x(k) - x^k)l + 1, k) f W(k) 

- W B (k) ] } 

l*(k+l, k)[x(k) - x«(k)] + T(k + 1, k)[W(k) - W B (k)]} T ) 

* (k + 1 , k) P(k) $ T (k+1, k) + *(k+ 1, k) K £ [7c(k) - 5c„(k)] 

[ W(k) - W B (k)l T ] r T (k + I, k) + r (k+ 1, k) E ( [W(k) - W.(k)l 

[ x(k) - x B (k)l T } 4> T (k + 1, k) + r<k+l, k) Q(k)T T (k+ 1, k) 

(B41) 

To evaluate the middle terms in the above, (B39) is subtracted from (B 20) 
to give 

x(k< 1) - x.(^ 1) $(k'l, k) fx(k) - x,(k)] -» 

' r<k *■ 1, k) ( W(k) - W«(k)] (B 42) 
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Let 


W(k) = W(k) - W,(k) 


x (k) = x (k) - 3?a(k) and 

Then 

x(k+l) = ♦(k + l, k) x(k) + T(k + 1, k) W(k) for k = 0, 1, ... 

Setting j = 0 in equation (B38), 

k 

x(k) = $(k, 0) x(0) + ^ $(k, i) T(i, i-l)W(i-l) (B43) 

i=l 


In analogy to (B 43) 

k 

x (k) = <fr (k, 0) x(0) <f>(k, i) T(i, i -1) W(i- 1) 

i=i 

and, therefore, 

E[x(k) W T (k)] = $(k, 0) E[x(0) W T (k)] + 

k 

+ X] i) r<I. i- 1) E[ W(i - 1) W T (k)] 

1=1 

W(k), k = 0, 1, ... is a zero mean Gaussian white sequence independent of the 
zero mean Gaussian random n vector x (0), i .e. , 

E [ x(k) W T (k)] - 0 

and 

E [ W(i - 1) W T (k)] = 0, i-1 i k 

Therefore, reverting to (B 41) 

P(k+ 1) - 4>(k + l, k) P(k) * T (k+l, k) + 

+ T (k+ 1, k) Q(k) r T (k+ 1, k) for k 0, 1, ... 

(B44) 

This equation gives the evolution of the covariance matrix for the discrete 
version of the stochastic process . 
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If t corresponds to time point k and t+At to k + 1 with At > 0, and 
considering that the analog of Q(k) for the continuous case needs to be 
obtained as Q(t)/At [Meditch, 1969], we have 

P(t+ At) = $(t+At, t) P(t) ® T (t+At, t) + 

+ r (t+At, t)|^ r T (t+At, t) 

Substituting for $ (t + At, t) and T (t + At, t) and expanding as in equations (b27) 
and (B 28), 

P(t+ At) = [I + A(t)At] P(t) [I+A(t)At] T + [G(t)At] ^p[G(t) At] T = 

= P(t) + A(t) P(t) At + P(t) A T (t)At 
+ G(t) Q(t) G T (t) At 

Transposing P(t) to the left and taking the limit as At ■+ 0 gives the matrix 
differential equation 

P = A(t)P + PA T (t) + G(t) Q(t) G T (t) for t * to (B45) 

with the initial condition P(to). This equation describes how uncertainty 
propagates in the system dynamics. The solution of this is of the form 

t 

P(t) - •(t, tb) P(tb) $ T (t, to) + f *(t, T) G(T) Q(T) G T (T) 4> T (t, T)dT 

to (B46) 

Equation (B46) is not of much use in the numerical computation of P(t) in 
general, as $(t, T) needs to be determined first. The numerical integra- 
tion of (B45) is preferred to evaluate P(t) . The last term in the right member 
of (B46) represents the contribution of state noise to P(t). 

7. Probabilistic Description of Measurements 
Recalling the observation state relation in the form 
y (t) H(t) x(t) + V(t) fort 5 to (B 47) 
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as in normal adjustment procedures, it is assumed that the measurement 
error T(t), t ^ to is p-dimensional white noise with E { V(t) ] = 0 and has 
covariance 

E {[V(t)][V(T)l T ] = R(t)6(t-r) for all t, T ^ to (B48) 

R(t) is continuous and positive definite. 

8. Classical Sequential Estimation 

With the mathematical details given above, the next step is to take up 
the problem of estimation. The models used in the experiment are special 
cases of the formulations studied above. In classical sequential estimation, 
noise is compensated for only as indicated in Section 4.1, and no effort is 
made to model the parameters of the state noise. This estimation procedure 
has been given here though not used specifically in this study since the pro- 
cedure for adaptive filtering follows from this procedure and clarifies the 
concepts underlying the Kalman filter algorithms. 

With G(t) taken as an identity matrix, the linear first-order differ- 
ential equation for the state deviation in (B6) becomes 

£(t) = A(t) x (t) + W(t) (B49) 

It is assumed that t s - to, the initial state X (t 0 ) ~ So. and the state noise W(t) 
is an n-vector with the properties 

E { W(t) } = 0, E { W(t) W(s) T ] = Q(t) 6(t- s) 

where Q(t) is the known n xn covariance matrix of state noise. X 0 

will not be perfectly known, and consequently the true solution X(t) will differ 

from the nominal solution X * (t) obtained with a specified initial state Xo. 

As a result, observations are taken to improve the estimate of X(t). 

The observation vector Y t , at discrete times ti, is related to the 
state by the functional form Yi - 0(Xj, ti) + Vi, i - 1, . . . , k - 1, k. 

Vi is the p-vector of observation errors with properties 
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E (Vt) = 0 


eCViVjM - Rifijj 


E {X(ti) V / j 


0 


Rt is a positive definite matrix which is assumed to be known. Lineai'ization 
around the computed observation with respect to the nominal state "X*(t), 
yields the linear form 

y,(t) = Hi x + Vi (B 50) 


where 


Hi 



* 

i 


which corresponds to equation (B47) given earlier. 

The solution of (B49) as seen from equation (B 3*1) is of the form 


x(t) = $(t, to) x(to) 


'/ 

to 


$(t, r) \V(T) dr 


(B 51) 


To commence the estimation, say, from a discrete time t k - 1 , the nominal state 

* JX 

X (t k _ I ) and estimate of deviation x (t*-*) are provided at time t*-!, along with 
an associated n x n conditional covariance matrix of the state 

A A A 

P*-a E { (Xk-x - Xk-x) (Xit-i ~ X*-i) T ! Y k -x ] (B52) 

where Y k -i implies conditioning on all observations from Y a through Y k -i. 

A A A 

With this understanding P k _i can be denoted by P(k-l/k- 1) and X*-! by 
/\ 

X(k - 1 A - 1). To commence the algorithm at k - 1, the initial estimates of 

A 

Po, even before the first observation vector is processed, is needed. In this 
study these estimates have been taken from the analyses in Chapter 4 . 

Similarly for k 1, the initial estimate of the state X 0 is taken from 
the broadcast ephemeris . 

Commencing at an arbitrary discrete time t*-x, a procedure is con- 
sidered for utilizing the observations Y k at the next discrete time t* to improve 
tlu* estimate of the state. This is done in the following steps: 
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Step 1 : 


Given the nominal state X * (t*-i) and the deviation x (tv-i) 

= F (X *, t), X*(t,c-i) = X(t*-i) 

and 

$(k, k-1) = A(k)$(k, k-1), $(k-l, k-1) = I 

are integrated numerically. The solutions will yield 5T *(tic) and 4>(k, k-1). 
The predicted state deviation 

£ (k/k - 1) = $(k, k-1) #(k-l/k-l) (B53) 

is then computed based on the concept of propagation of the mean, given 
in equation (B39). 

Step 2 : 

The predicted error covariance matrix of the state P(k/k- 1) is com- 
puted from 

P(k/k- 1) = 4>(k/k- 1) P(k-l/k-l) $ T (k/k- 1) + Q(k) 


where 


Q(k) = 


f ^(tic, T) Q(T) <& T (tk, T)dr 
1 vide equation (B 46) 


Analytical expressions for this have been derived in [Myers, 1973]. 


Step 3 : 

The observation vector Y(t k ) is now processed for estimation by first 
computing the observational deviation 

* V>- 0 ( *(M +i(k/k - ,H (B54) 

and the associated H matrix for the observation state relation 
(tk) 


= 'W‘ k/k) * v «w 


(B55) 
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in which x(kA) and V.. . are the unknown state deviation and the observational 

(Me) _ _ 

error. Let H be denoted as H , , and v., as y„ 

0*) (k) • (W) (k). 

Under the assumption that the estimate is a linear function of observa- 
tions, and using the principle of minimum variance of the estimate as in 
normal least squares procedures, the following expression for the estimate of 
state at step k, after processing observations at t*, is obtained 

A A A 

x (k/k) = x(k/k-l) + K(k) [y - H (R) x (k - 1/k - 1)] (B56) 


where K(k) is the n x p Kalman gain matrix given by 


K (k) J>(k/k " l) H (k) fH (k) P (k/k - 1) H (k) + R (k) 1 


(B57) 


This matrix is also called the filter pain. In (B 57) the matrix to be inverted is 
of size p x p, the number of observations at t*. The observation cova:- ianee 
matrix R„ . is assumed to be positive definite to ensure that the matrix in the 

(k) 

brackets can be inverted. The new estimate of the state at tv is 


X 


- x (U + *** 


Step 4 : 


The new error covariance of the state based on the observations up to 

and including Y .. . will be 
(*) 

P(k/k) - [I - K H 1 f>(k/k-l) (B58) 


and the algorithm is repeated for the next observation at tv+i. To minimize 
the linearization errors, the current best estimate of the state is reinitialized 
as the nominal state vector for the next step. This procedure makes 



x * (W ♦ j(k/k, 


the reinitialized nominal state for the next step, setting state deviation in 
Step 1 above to zero. This also simplifies computations. 
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It can be seen that the above estimation procedure which is at times 
known as first-order, nonadaptive filtering procedure with state noise com- 
pensation is very similar to a normal sequential least squares procedure, as 
applied to a dynamic situation. The inclusion of state noise is comparable to 
the procedure of weighting parameters in a solution. For example, we can 
compare equation (B 56) above to equation (27) of [Uotila, 1975] . This pro- 
cedure treats state noise as a white noise process and takes cognizance of it 
in the adjustment procedure by assuming its statistical properties. 
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